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ABSTRACT 

The formalism to compute the geometrical and topological one-point statistics of 
mildly non-Gaussian 2D and 3D cosmological fields is developed. Leveraging the 
isotropy of the target statistics, the Gram-Charlier expansion is reformulated with 
rotation invariant variables. This formulation allows to track the geometrical statis- 
tics of the cosmic field to all orders. It then allows us to connect the one point statistics 
of the critical sets to the growth factor through perturbation theory, which predicts 
the redshift evolution of higher order cumulants. In particular, the cosmic non-linear 
evolution of the skeleton's length, together with the statistics of extrema and Euler 
characteristic are investigated in turn. In 2D, the corresponding differential densities 
are analytic as a function of the excursion set threshold and the shape parameter. 
In 3D, the Euler characteristics and the field isosurface area are also analytic to all 
orders in the expansion. Numerical integrations are performed and simple fits are 
provided whenever closed form expressions are not available. These statistics are com- 
pared to estimates from N-body simulations and are shown to match well the cosmic 
evolution up to root mean square of the density field of ~ 0.2. In 3D, gravitational 
perturbation theory is implemented to predict the cosmic evolution of all the rele- 
vant Gram-Charlier coefficients for universes with scale invariant matter distribution. 
The one point statistics of critical sets could be used to constrain primordial non- 
Gaussianities and the dark energy equation of state on upcoming cosmic surveys; this 
is illustrated on idealized experiments. 



1 INTRODUCTION 

Random fields are ubiquitous in physics. In cosmology, the large scale distribution of matter (LSS) and the sky-maps of 
the polarized Cosmic Microwave Background (CMB) radiation are described as, respectively, 3D and 2D random fields. In 
primordial theories where initial seeds for cosmic structures have quantum origin, the fields of initial density fluctuations are 
close to Gaussian. Subsequent evolution develop further non-Gaussian signatures via non-linear gravitational dynamics. Hence 
these signatures might be used to constrain the dark energy equation of state via 3D galactic surveys, or shed light on the 
physics of the early Universe. The non-Gaussianities can be accessed through the higher order moments of the field, which are 
generally difficult to estimate directly in real-life observations, due to their sensitivity to very rare events. The geometrical 
analysis of the critical sets of the field provides more robust measures of non-Gaussianity and has become an active field of 
investigation (Gott et al. 2007; Park et al. 2005). Another active area of research is centered around the use of the bi- or 
trispectrum for evaluation of the non-Gaussian effects (Scoccimarro et al. 1998; Scoccimarro, Sefusatti & Zaldarriaga 2004; 
Sefusatti & Komatsu 2007; Liguori et al. 2010; Sefusatti, Crocce & Desjacques 2010; Nishimichi et al. 2010). This paper 
presents statistical tools that may help to detect unique signatures of the fundamental processes in our Universe. 

Several statistics have been used to characterize the morphology of the density field of the large-scale structures. Peaks 
are one of the first components to have been studied in detail (Bardeen et al. 1986) . Maxima of the density fields are indeed 
considered as the location where gravitational collapse occurs to form the halos where galaxies form. Their distribution is 
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thus a key ingredient in semi-analytical models and for predictions of the so-called halo model (Mo & White 1996; Cooray & 
Sheth 2002). In 2D, the statistics of peaks of the cosmic microwave background has been used to constrain its Gaussianity 
(Bond & Efstathiou 1987; Larson & Wandelt 2004; Hou et al. 2010b; Pogosyan, Pichon & Gay 2011). 

The Euler characteristic and the genus, which only differ by their normalization, have also been thoroughly studied, 
since they describe in part the topology of the field. They provide a different approach to the usual geometrical probes 
such as the power-spectrum. The genus has been used in redshift surveys (Gott, Dickinson & Melott 1986; Vogeley et al. 
1994; Canavezes et al. 1998; Hikage et al. 2003; Park et al. 2005). It has also been applied to the CMB to constrain the 
Gaussianity of primordial fluctuations, with results from COBE (Smoot et al. 1994; Colley et al. 1996; Kogut et al. 1996) and 
WMAP (Colley & Gott 2003; Gott et al. 2007). The genus can be defined in a probabilistic framework where its theory is 
more easily derived. Gaussian predictions have been known for a long time (Doroshkevich 1970; Adler 1981) and have been 
extended by Matsubara to non-Gaussian fields arising from gravitational non-linear dynamics (Matsubara 1994) and to proper 
accounting of redshift distortions (Matsubara & Suto 1995; Matsubara 1996). More generally, the Euler characteristic is one 
of the Minkowski functionals (a complete set of morphological descriptors) (Mecke, Buchert & Wagner 1994a) and have been 
combined to describe the shape of structures (Sahni, Sathyaprakash & Shandarin 1998). This paper is an extension of this 
approach to other characteristics and higher orders. 

The skeleton has been introduced more recently (Novikov, Colombi & Dore 2003). It traces the filamentary structure of a 
field, extracting the "ridges" linking the maxima through saddle points. It has been applied in 2D for the CMB (Eriksen et al. 
2004; Hou et al. 2010a) and in 3D for the study of the large-scale structure (Sousbie et al. 2008; Aragon-Calvo et al. 2007; 
Platen, van de Weygaert & Jones 2007). Most of the efforts have focused up to now on improving the algorithm (Sousbie et al. 
2007; Sousbie, Colombi & Pichon 2009; Sousbie, Pichon & Kawahara 2011) and our theoretical knowledge: Pogosyan et al. 
(2009) gives a complete theoretical description of the properties of the skeleton for Gaussian fields. One of the motivations of 
the present paper is to extend this description of the skeleton beyond this specific Gaussian case. Indeed, mildly non-Gaussian 
random fields play an important role in cosmology: they can be present in the primordial fluctuations and open a window 
to the physics during the inflation phase (Bartolo et al. 2004), and they also arise e.g. from the non-linear dynamics of 
gravitational clustering (Bernardeau et al. 2002). 

In this paper we develop the formalism for computing the geometrical and topological one-point statistics for mildly 
non-Gaussian 2D and 3D cosmological fields. Our formalism is based on the Gram-Charlier expansion of the Joint Probability 
Distribution of the field and its derivatives using rotation invariant variables which allows to analyse the isotropic one-point 
statistics to all orders in cumulant expansion. In the following we will predict the evolution of quantities such as the Euler 
characteristic, the density of extrema and the differential length of the skeleton as a function of root mean square of the cosmic 
density field, a. The evolution of these geometrical critical sets can thus be linked to the cosmic evolution of the underlying 
field. The next generation of redshift surveys, such as Big-BOSS or EUCLID and WFIRST, 1 aim at constraining the equation 
of state of dark energy by mapping millions of galaxies. For example, Zunckel, Gott & Lunnan (2011) showed how using the 
genus in Big-BOSS could give a precision of 5% on a constant equation of state. Similarly, Yang et al. (2011) argue that a joint 
powerspectrum and extrema counts probe has approximately twice the cosmological sensitivity of powerspectrum estimation. 
As an alternative route, we will show below, for instance, that the length of the skeleton per unit volume within some range 
of some density threshold can be written as L(z) = Z/ ' + o(z)8L^ where £A ' is the Gaussian prediction and SlS 1 ' is the 
non-Gaussian prediction at first-order in a. By knowing these two functions and measuring L(z), <j(z) can be extracted and 
the skeleton could thus be used to measure the growth factor of the fluctuations in a manner which might be less sensitive 
to biases. Indeed, a common feature of the critical sets that we will build below is that they are formally independent of any 
monotonic re-mapping of the underlying field. 

The structure of the paper is the following. In Section 2, we introduce the Gram-Charlier expansion and the rotational 
invariants used to describe an isotropic field. We obtain the expression of the expansion in two and three dimensions. In 
Section 3, these expressions are used to compute a sequence of topological invariants in 2 and 3D: the Euler characteristic, 
the distribution of extrema and the length of the skeleton. The probability distribution of the field and the statistics of its 
iso-contours are also given. Section 4 presents applications to CMB non-Gaussianities and to the cosmic evolution of the 
large scale structures. In particular, Section 4.2.2 shows how to predict the cumulants that appear in these expansion using 
perturbation theory when cosmic gravitational clustering is responsible for their growth, Section 4.2.3 compares the results to 
measurements on numerical realizations and Section 4.2.4 describes an idealized dark energy probe experiment. Appendix D 
details the link between the order in the Gram-Charlier expansion and that corresponding to gravitational perturbation theory, 
introducing the so-scalled Edgeworth expansion. It is then applied to the computation of the <r 2 -order correction to the Euler 
characteristic in Appendix E. The results from perturbation theory to scale invariant power spectra are given in Appendix F. 
Appendix G presents the algorithm implemented to measure the Euler characteristic. 



http://bigboss.lbl.gov/, http://sci.esa.int/euclid/, http://wfirst.gsfc.nasa.gov/ 
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2 ASYMPTOTIC EXPANSION OF THE JOINT PDF OF THE FIELD AND ITS DERIVATIVES 

The expectation of topological and geometrical statistics like the Euler characteristic, extrema counts, the differential length 
of skeleton requires the knowledge of the joint probability distribution (JPDF) of the field and its derivatives up to second 
order (Adler f981). In this paper we advance the formalism that considers these statistics in non-Gaussian regime treated as 
the correction to the Gaussian limit. We start with developing the relevant version of non-Gaussian JPDF of the field and its 
derivatives, following closely the exposition in Pogosyan, Gay & Pichon (2009). Starting with some notations, for a given field 
p, we define the moments a 2 = (p 2 ), o\ = ((Vp) 2 ), cr| — ((A/ 9 ) 2 ) - Combining these moments, we can build two characteristic 
lengths, Ro = a/ai and i?» = tri/az, and the spectral parameter 7 = o\l(ooi). We choose to normalize the field, p and its 
derivatives to have unit variances: x — p/a, Xi — l/<7iV;p, and Xij = l/craVjV^p. 



2.1 The Gram-Charlier expansion of the JPDF 

A statistically homogeneous N dimensional random field x is fully characterized by the joint one-point distribution function 
(JPDF) of its value and its derivatives -P(x), x = (x,Xi,Xij,Xijk, ■ ■ •)• The simplest case, frequently arising because of the 
central limit theorem (e.g. simple inflation models), is the Gaussian JPDF: 

G(x) = (2tt)- n/2 \C\- 1/2 exp(--x ■ CT 1 ■ x) , 

2 

where C = (x T ■ x) is the covariance matrix of the x variables and we have considered that all the variables are defined 
as having zero mean and unit variance (Cm = 1 for each i). For a general non-Gaussian field, the JPDF can be expanded 
around a Gaussian arranged to match the mean and the covariance of the field variables x, using a multivariate Gram-Charlier 
expansion (Chambers i967; Juszkiewicz et al. 1995; Amendola 1996; Blinnikov & Moessner 1998): 



P(x) = G(x) 



00 

1 + Ei Tr Kx n >GC-h„(x)] 



(1) 



The correction to the Gaussian approximation is the series in Hermite tensors 

h n ( x ) = (^i)"G~ 1 (x)a n G(x)/ax™ , 

of rank n. These tensors are the orthogonal polynomials related to the Gaussian kernel G(x), which allows to construct the 
so called Gram Charlier coefficients from the moments: 

(x") G0 = (h„(x)>. (2) 

The statistics of the Euler characteristic (or equivalently the genus) and extremal points (Adler 1981), and the basic description 
of the critical lines of the field in the stiff approximation (Pogosyan et al. 2009), require to know the JPDF only up to the 
second derivatives. In the Gaussian limit, in this case, the only non-trivial covariance parameter is the cross-correlation 
between the field and the trace of the Hessian 7 = — {xTr(xij)) (Bardeen et al. 1986; Pogosyan et al. 2009). 



2.2 Rotational invariants of the field and its derivatives 

The JPDF in the form of equation (1) is not ideal to study critical sets statistics since the coordinate representation does 
not take advantage of the isotropic nature of the statistical descriptors. The pioneering works of Matsubara (1994, 2003), 
where the first correction to Euler characteristic was computed, demonstrate the arising complexities. Instead we developed 
the equivalent of the Gram-Charlier expansion for the JPDF of the field variables that are invariant under coordinate rotation 
(Pogosyan, Gay & Pichon 2009; Pogosyan, Pichon & Gay 2011). Such distribution can be computed via explicit integration 
of the series (1) over rotations. However we obtain it here directly from general principles: the moment expansion of the 
non-Gaussian JPDF corresponds to the expansion in the set of polynomials which are orthogonal with respect to the weight 
provided by the JPDF in the Gaussian limit. Thus, the problem is reduced to finding such polynomials for a suitable set of 
invariant variables. 

The rotational invariants that are present in the problem are: the field value x itself, the square modulus of its gradient, 
q 2 — Yli x i and the invariants of the Hessian, the matrix of the second derivatives x%j. A rank N symmetric matrix has 
N invariants with respect to rotations. The eigenvalues Xi provide one such representation of invariants, however they are 
complex algebraic functions of the matrix components. An alternative representation is a set of invariants that are polynomial 
in Xij, with one independent invariant polynomial per order, from one to TV. A familiar example is the set of coefficients, 
7 S , of the characteristic equation for the eigenvalues, where the linear invariant is the trace, Ji = Ai, the quadratic one 
is I2 = ~}2 i< j ^i^j an d the AT-th order invariant is the determinant of the matrix. In = FL-^- Aiming at simplifying the 
JPDF in the Gaussian limit 2 (e.g. Doroshkevich 1970; Pichon & Bernardeau 1999; Pogosyan et al. 2009) we use their linear 



2 By Gaussian JPDF we mean the limit of the Gaussian field x. Note that even in this limit, the rotation invariant variables themselves 
are not Gaussian-distributed, if they are non-linear in field variables x. 
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combinations, J s , 

Ji = h , Js>2 = I{ - ^2 7Z _ i i~ Pi p > ( 3 ) 

p=2 1 

where J s >2 are (renormalized) coefficients of the characteristic equation of the traceless part of the Hessian and are independent 
in the Gaussian limit on the trace Ji. 

Our choice of normalization of the field and its derivatives gives simple relations for the linear and quadratic in the field 
moments of the invariant variables: 

(0 = o, (Jr) = o, (CJi) = o, (c 2 ) = l, <Ji 2 ) = l, ( q 2 ) = 1, (J 2 ) = 1 , (4) 

i.e., the linear invariant variables have zero means and unit variances and the quadratic ones have unit means. Since we expand 
JPDF around the Gaussian form that is arranged to match all these terms, the relations (4) are retained in non-Gaussian 
case. 



2.3 The Gram Charlier non-Gaussian joint probability distribution in two dimensions 



£11X22 — %2 ■ 



The rotation invariant formalism is straightforwardly adapted to 2D random fields, which are of interest for e.g. CMB and 
weak lensing studies. There are just two invariants that describes the Hessian in 2D 

h = Ti(xn) = m + X22 , and I2 = det \xa\ = 

which in the Hessian eigenframe have the form 

h = Tv(xij) — Ai + A2 , and I 2 = det \xt 

We shall also use the orthogonal set 



Ai A2 



(5) 
(6) 
(7) 

Introducing £ = (a; + 7Ji)/\/l — -y 2 in place of the field value x we can build the 2D Gaussian JPDF G 2 d(C,9 2 , Ji, J a ), 
normalized over d£dg 2 dJidJ2, which has a fully factorized form: 



Ji — h, and J2 = J 2 — 4/2 = (xii — x 22 ) 2 + 4a; 2 2 . 



G 2 d(C, q 2 , Ji,J2)dC,dq 2 dJ 1 dJ 2 = ~- exp 

Z7T 



1 A 2 2 

'2 C q 



1 7 2 / 
2 Jl ~ 32 



dC^dq dJidJ'2- 



(8) 



From this Gaussian kernel, we can build the basis of orthogonal polynomials associated with it. Since the variables are 
independent, the polynomials will not be coupled and we can consider each variable separately. For variables with a normal 
distribution (£ and Ji), the orthogonal polynomials are the "probabilists" Hermite polynomials, which follow the orthogonality 
relation: 



/2n 



exp 



H„(x)H m (x)dx = n\ 5™ 



For the quadratic variables (q 2 and J2), the corresponding polynomials are the Laguerre polynomials: 



r 

Jo 



exp(—x)L n (x)Lm(x)dx = 8^ 



(9) 



(10) 



The non-Gaussian rotation invariant JPDF can then be written in the form of the direct series in the products of these 
polynomials: 3 



J , 2d(C,9 ,Ji,Ja) = G ! 2D 



00 i+2j+k+2l—n 
n— 3 i,j,k,l 



M) J+i 

i\ jl fe! l\ 



Cq 2j Jl k J2 l ) Hi (C) Lj (q 2 ) H k (JxjL, (J 2 



(11) 



where ^2l +2 j!'l k+21 ~™ stands for summation over all combinations of non-negative i, j,k,l such that i + 2j + k + 21 adds to the 
order of the expansion n, defined as the number of fields. This arrangement sorts the series in the order of increasing power of 
the field variables. 4 When counting the order of a term, j and I enter multiplied by two, since each of q 2 and J2 contributes 
a quadratic in the field quantity. 

The orthogonality of the polynomials (Abramowitz & Stegun 1970) allows us to extract the coefficients by projection: 



, . r-\-OG z' + oo /--t-oo /- + 00 

Hi (C) Lj (q 2 ) H k (Ji) U {Ja)) = dC / dg 2 / dj 1 / dJ 2 P 2 D(C, q 2 , Ji, J 2 )^ «) Lj {q 2 ) H k (Ji) U 

1 J -00 JO J -00 Jo 



(J 2 



(~l) J+i 



(12) 



Cq 23 Ji k J2 



3 The orthogonality expressions and the choice of normalization for the Gram-Charlier coefficient (equation (12)) explain the quadruple 
factorial at the denominator together with the (— iy +l factor. 

4 This ordering should not be confused with the order given by perturbation theory. See Appendix D for further explanations. 
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We can thus express the Gram-Charlier coefficients as combinations of moments: 

(CV Vf4) oc = J~^f+i { H i (0 U (l 2 ) H k (Ji) U (J 2 )) . (13) 

The chosen convention for the prefactor ensures that when the polynomial expressions are expanded, the highest power term 
is given directly by the corresponding moment: 

(cv i #4) oo =(cvvi4)+... . (14) 

For example, 

<C 3 ) GC = (# 3 (C)> = <C 3 -3C> = <C 3 >-3<C> = vC 3x 



(q 2 Jl) GC = (L 1 (q 2 )H 1 (J 1 )) = -{(l-q 2 )J 1 ) = { q 2 J 1 )-{J 1 ) = (q 2 J 1 ) 



(J 2 J 2 ) oc = ( J ff 2 (Ji)Li(J 2 )) = - ((J 2 - 1)(1 - J 2 )> = <J?J 2 ) - (J?) - (Ja> + 1 = (Ji 2 J 2 ) 
where we used our normalization convention equation (4). Consequently, the n — 3 Gram-Charlier coefficients simply coincide 
with the third-order moments, but the higher order ones are given by the combination of moments. 

The combinations of moments that constitute the Gram-Charlier coefficients vanish in the Gaussian limit. Thus they can 
be expressed via the cumulants of the field that have the same property. For the variables that are linear in the field, e.g. 
C, Ji, the ordinary cumulants arise. However, the variables q 2 and J 2 are not Gaussian distributed in the limit of the Gaussian 
field due to their quadratic nature, and their own cumulants do not vanish in this limit, e.g., (j|) cllm = (jf) — (J 2 ) 2 — > 1. To 
take this effect into account we introduce the notion of "field cumulants". A "field cumulant" of a polynomial variable is the 
cumulant of the expression in field variables, x, Xi,Xij, that constitutes it. By definition it vanishes in the limit of the Gaussian 
field. For example, (jf )f c = ( [O^n — X22) 2 + 4x 2 2 ] 2 )c = {Ji} — 2( J 2 ) 2 — ¥ 0. In the rest of the paper we will always mean "field 
cumulants" and will use the subscript "c" to tag them. Expressing the Gram-Charlier coefficients via the field cumulants, we 
find that for the orders n — 3, 4, 5 they coincide, e.g. 

(CV = " (H2(t)Li(q 2 )Hi(Jl)) = - <(C 2 - 1)(V + 1) Jl> = «V Ji) - (t 2 Ji) - (q'Ji) = (CV Ji) c ■ (15) 
This property follows from the link between the Gram-Charlier and the Edgeworth expansions (see Appendix D). For the next 
orders (n > 6), the Gram-Charlier coefficients are more complex combinations of cumulants, as explained in Appendix E. 

2.4 The Gram Charlier non-Gaussian joint probability distribution in three dimensions 

In 3D, the polynomial invariants of the matrix Xij, as written in the coordinate frame, are 

h = Tr(xij) = xu +X22 +X33 , h = X-1XX22 + X22XZZ + x^xyj, — x\ 2 — x\z — x \i , (16) 
I 3 = det \xij\ = x-i\X22Xzz + 2x\2X2zXya — xnxl 3 — X22x\ z — x 33 a; 2 2 . (17) 

The same invariants expressed through eigenvalues are 

h = Tr(asii) = Ai + A 2 + A 3 , h = A1A2 + A2A3 + A1A3 , and J 3 = det \x tJ \ = AiA 2 A 3 . (18) 

This variables can be made partially independent using the linear combinations 

q 07 

Ji=/i, J 2 = 7i 2 -3/ 2 and J 3 = i? - -hh + y h. (19) 
In iV dimensions, the Gaussian JPDF, Gnd(C q 2 , Jit Js>2), retains complete factorization with respect to £, q 2 , Ji 



iv^*„JV-2 r t 2 N 2 T 2 



(%)*q»-< r ,2 , r „ 2 T2 

Gnd= 2.r[f] exp 



so in the moment expansion this sector always gives rise to Hermite Hi((), Hh(Ji) and generalized Laguerre L\ N ~ 2)/2 {Nq 2 /2) 
polynomials. However, the distributions of the rest J s > 2 contained in Q{J S >2) are coupled. Specifically, in 3D, Q{J2, Jz) 
exp[— I J 2 ], and J 3 is distributed uniformly between —j'^ 2 and j\ 



3/2. 



Gzr>((,q 2 ,Ji,J2,Jz) 



8tt 2 

with the ranges that variables span summarized in the normalization integral 

,3/2 



1 > 2 3 2 1 7 2 5 
2 C 2 q 2 Jl " 2 J2 



(21) 



dC / dq 2 / dJi / dJ 2 / dJzGzD((,q 2 ,Ji,J2: 
-00 Jo i-00 Jo J-jI' 2 



Js) = 1 • (22) 



Let us denote the orthogonal polynomials in the variables J 2 and J 3 , Fl m (J 2 , J 3 ), where Z is the power of J 2 and m is the 
power of J 3 . They obey the orthonormality condition 

.3/2 

dJ 2 / ; dJ 3 ^Rm(J 2 , Jz)F Vm ,(J 2 ,J 3 ) = <5| C • (23) 



3/2 
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We shall not give the full theory of these polynomials here, but note that one can construct them by a Gram-Schmidt orthog- 
onalization procedure (say using the Mathematica Orthogonalize routine using the dot product defined by equation (23) 5 ) to 

any given order. Two special cases Fio = \J (3+2i)l M^ 3 ^ (I an( ^ -^ 01 ~ vli^ 3 are sumc i en t to obtain the general expression 
for the Euler characteristic of the excursion sets of the field to arbitrary order, and to calculate the critical point and skeleton 
statistics to quartic order. Hence, we write the moment expansion for invariant non-Gaussian JPDF, P3d(C </> J\,Ji, J3) as a 
series in the power order of the field, n in the form 



G 3 



n— 3 i,j,k,l 



00 i J r2j + k J r 3=n 

E E 

n — 3 i,j,k 

00 i+2j + k + 2l+3m=n 

+ E E 

n= 5 i,j,k,l,m — l 



i! (1 + 2/)!! k\ (3 + 20!! 
(-1) J 3 3 x 25 



Ci 23 7 k j I 
q Ji Ji 



.(1/2) 
"J 



i\ (1 + 2j)l\ k\ x 21 



H> (C) L 

Hk (Jl) J3 



3 2 
2 q 



(3/2) 



5 1 

2 J2 



+ 



3 2 
7,1 



i\ (l + 2j)!! fc! 



z-i 23 -r k t I r Tti 

C q Ji J 2 J3 



^(C)4 1/2) 



3 2 



H fe (J 1 )F, m (J 2 ,J 3 ) 



(24) 



The first non-Gaussian correction term contains J2 but not J3, the second is linear in J3 and contains no J2, while the 
last one contains all the remaining combinations of J2 and J3. In this last term we leave the normalization coefficients c; m 
undetermined, but coi = 0, so the first non-zero contribution, I = m = 1, is of the fifth power of the field and proportional to 
J2 J3. As in 2D, the orthogonality can be used to express the coefficients in terms of moments. For m — 0: 



z-i 23 t k t I 



and for 1 = and m = 1: 



+i 2J t k t 

C Q Ji J3 



(3/2)' (5/2) ; 
_ (-1)^1 



(1/2) 



Hi (C) i 



ffi(04 1/a) 



3 2 
2^ 



(3/2) 



3 2 



H k {Ji)L\ 



Hk (Ji) J3 



:J 5 



(25) 



(26) 



(3/2)' 

Equations (11) and (24) are the main results of this section. They fully characterize the joint probability function of a field 
and its derivatives in two and three dimensions, up to second order derivatives, to an arbitrary order in the moment expansion 
(see also Pogosyan, Gay & Pichon (2009); Pogosyan, Pichon & Gay (2011)). 



3 NON-GAUSSIAN CRITICAL SETS IN 2 AND 3 D 

Given formulas (11) and (24) one can easily compute any rotation-invariant statistics that depend exclusively on these 
descriptors of the field. In particular, one can compute the filling factor of the field excursion set (Section 3.1), the length of 
field isocontours in 2D and the area of isosurfaces in 3D (Sections 3.2.1 and 3.3.1), the Euler characteristic of the excursion 
sets of the field (Sections 3.2.2 and 3.3.2), the density of extrema (Sections 3.2.3 and 3.3.3), and the properties of the critical 
lines that describe the skeleton of the field (Sections 3.2.4 and 3.3.4). 



3.1 The filling factor of high excursion set 



The filling factor, which is the volume fraction occupied by the field with values exceeding the threshold v, presents the 
simplest case when one considers only the field and not its derivatives. The JPDF, equations (11) and (24), is then reduced 
both in 2D and 3D to the same expression: 



P(x) = G(x) 



i+E 



[ H n (x) 



(27) 



This is the classical one-variable expression for the Gram-Charlier expansion. This expansion does not have good convergence 
properties (Blinnikov & Moessner 1998), but it is very simple to implement. When the origin of non-Gaussianities is known and 
can be used to predict the amplitude of each coefficient, this expansion can be re-summed to a more physically motivated one. 
For example, as shown in Appendix D, when the non-Gaussianities give rise to a hierarchy of cumulants such that {x"} c oc a p ~ 2 
(notably in the case of gravitational evolution), the re-summation to a series in terms of a leads to the Edgeworth expansion 
(Juszkiewicz et al. 1995): 

2 



P(x) = G{x) 



1 + 



H s (x) + 



24 



HJx) + 



72 



H 6 (x] 



+ ■ 



(28) 



for instance, the orthogonal polynomial of order one in J2 and J3 is proportional to (5 J2 — 11) J3. 
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From the non-Gaussian PDF given by equation (27), the filling factor, f(u) immediately follows 

fW) = P(x)dx = ^Erfc + G{v) f] i^L Hn ^(u). (29) 

3.2 Non-Gaussian critical sets in two dimensions 

3.2.1 The length of 2D isocontours 

The average length of the isocontour x = v is one of the two Minkowski functionals (Schmalzing & Gorski 1998; Mecke, 
Buchert & Wagner 1994b; Schmalzing & Buchert 1997) that characterize the level sets in 2D. In our formalism this is the 
simplest statistics for which the invariant variables appear non-trivially. To compute the average length per unit volume C of 
a ^-isocontour, 

£(!/) = —/ dx dq 2 P 2D (x,q 2 )5 D {x - v)q = — / dq 2 P 2D (v, q 2 )q , (30) 

-KO y_oo Jo H J 

one needs JPDF of the field and the magnitude of its gradient. From equation (11) 

1 2 / °° i + 2 J= n / -. \j \ 



yielding 



C( V ) = -^ e -£ ( T e^qdj + j^^" ^Uq V ) [°° e^ qL^dq 2 ) . (32) 

The identity J °° e -9 qLj(q 2 )dq 2 = —T(J — |)/(4T(j + 1)) allows us to write the length of the v isocontour to an arbitrary 
order in non-Gaussianity as 

£(,) = *e-* (l + 4= £ IT ( ~;r ?| (*V ) %)V (33) 

W 2^2% V 2 ^^a « ^ ! r + 1 ) X /gc 7 

This result is a generalization of the published results to first (Matsubara 2003) and second (Matsubara 2010) orders. Ap- 
pendix E2 reorders the expression (33) as an Edgeworth expansion. 

3.2.2 The Euler characteristic of non-Gaussian 2D fields 

The Euler characteristic of a manifold is defined as the alternating sum of the Betti numbers: x — bo — bi + 62 — . . . . The first 
Betti numbers have a simple interpretation: 60 is the number of connected components, 61 is the number of two-dimensional 
holes and 62 is the number of three-dimensional voids. It is also one of the Minkowski functionals (Mecke, Buchert & Wagner 
1994a) and in particular it is additive. For example the Euler characteristic is 1 for a ball and 2 for a sphere (which is a 
region surrounding a void). The Euler characteristic is an important quantity in topology and plays a central role in many 
different theories. For example, it can be linked to the geometry of the surface via the Gauss-Bonnet theorem or the extrema 
of a function in the Morse theory (Jost 2008). 

In the specific case of an orientable and closed two-dimensional surface, it can be shown to be topologically equivalent to 
a sphere with g handles. One can thus use g, called the genus, as an equivalent to the Euler characteristic. There is a direct 
relation between the Euler characteristic of a closed surface and its genus: x = 2 — 2g. 

In cosmology, the genus has been used to characterize the topology of the isocontours. However the definition generally 
used by cosmologists is that the genus is the number of holes minus the number of isolated regions, which leads to an offset 
of 1: for an isolated region, the genus is then defined to be the number of holes minus 1 (see e.g. the warning in Gott et al. 
(2008)). The advantage of this shifted definition is due to the fact that the relation between the Euler characteristic x an< i the 
"cosmological" genus g c in this definition is simply: x — ™2g c . One can therefore keep the intuitive definition of the genus (i.e. 
the number of holes) while inheriting all the good properties from the Euler characteristic (i.e. the Gauss-Bonnet theorem and 
the results from Morse theory). Table 3.2.2 gives the values of the genus (in the two definitions) and the Euler characteristic 
for basic sets. In this paper, we will focus on the Euler characteristic of the excursion sets, i.e. the volumes above a given 
threshold. This approach is equivalent but nonetheless different from looking at the genus of the isocontour surface. However, 
the difference only affects the prefactor since, for a 3D field, the Euler characteristic of the isocontours is simply twice the 
Euler characteristic of the corresponding excursion sets (Mecke, Buchert & Wagner 1994a): x(9M) = x(^) [l + (~ 
where d — dim(M). To compute the Euler characteristic of a random field, we use a central result from Morse theory (Jost 
2008). Morse theory makes the link between the topology of a manifold and the critical points of the functions defined upon 
it. It shows that the Euler characteristic is simply equal to the alternating sum of the number of critical points of index i (i.e. 
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Table 1. Values of the Euler characteristic and the different definitions of the genus, for basic topologies. 



where the Hessian has i negative eigenvalues): \ = l)* n »- F° r a two-dimensional random field, the Euler characteristic 

x{v) of the excursion set x > v is given by: 



(34) 



where n max (^), n BW x[y) and n m i n (^) are the number of maxima, saddle points and minima above the excursion threshold, v. 
The density of Euler characteristic can thus be computed by integrating the determinant of the Hessian, I2: 

1 

M 



dJi 



dC / dJ 2 P cxt (C, Ji, Ji)h, 



(35) 



where P ex t is the distribution function under the condition of zero gradient: 



Pext((,Jl,J2) = 



dq 2 P 2 D(C,<? 2 , Jl, J2)8d{xi)5d{X2)- 



(36) 



As demonstrated in Appendix Al, this integration yields: 



-Erfc 



V2j 



x(-oo) + 



1 



+ E E 

n — 3 i,j,k 



(-iy +k 



exp 



GC 



7f/iM + 

min(2,S;) 



•7 ) 



00 i+2j + k+2 = n 

+ E E 

n— 3 i>j>k 



E 

j+fc+1 



2! 7 



fe+2-2s 



7 Hi 



(-1) 

i! j! fc! 



s!(2 - s)l(k - s)\ 

\ / GC 



(37) 



Note that in the rare event limit, v > 1 or v < 1, equation (37) represents resp. the differential number of maxima or minima, 
as the contribution to equation (34) from other extremal points becomes negligible (Aldous 1989). In this limit, equation (37) 
yields both quantities to all orders. 

A comment is in order on the value of the Euler characteristic density evaluated at infinitely low threshold, x(~ °°)- 
In a finite, topologically non-trivial space 6 x(~°°) evaluates to the Euler characteristic of the manifold itself, divided by 
its volume. While the result of Morse theory in equation (34) is exact, the expression via the local estimate of the extrema 
densities, equation (35), does not contain this effect as {I2) = 0. The expansion in equation (37) shows the appearance of the 
boundary term ^Erfc(z/)x(— 00) when the moments of the form (Jlq 21 ) — (J 2 g 2! ) do no cancel exactly. We conjecture that 
such boundary term accounts for the topology of the manifold on which the field is considered. It might perhaps be one of 
the most important applications for the genus statistics of masked data sets with the mask introducing non-trivial topology. 

When the non-Gaussianities give a hierarchy of cumulants such that a cumulant of n fields is of order a n ~ 2 (notably 
for gravitational evolution), the Gram-Charlier series can be reordered as the series in powers of a (called an Edgeworth 
expansion, see Appendix D): 



x(y) 



-Erfc 
2 



x(-oo) + 



47TV27TP 



exp 



\) 7 2 #i(f) + E 



(38) 



6 The only topologically non-trivial space in 2D that is also maximally symmetric, i.e., homogeneous and isotropic, is the two-sphere 
S 2 , the case relevant to CMB studies. Due to the curvature of S 2 , the invariant second derivatives are expressed via covariant derivative 
tensor with mixed components, x-^ , J\ = x-^'' and J2 = J\ — 4|a3^ u '|. In all spaces where the global isotropy is broken, as, for example, 
in all topologically compactified fiat spaces, like torus, the invariant formalism gives an approximation which is more accurate smaller 
the R t scale is relative to the size of the manifold. 
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Figure 1. Left: Euler characteristic in 2D slices of a mildly non-Gaussian field (gravitational clustering, cr = O.f). Dashed line: Gaussian 
approximation, solid line: first order prediction using measured cumulants (following the prescription described in Section 4.2.3). The 
narrow shaded band corresponds to the 3-sigma dispersion for the measurement of the Euler characteristic in the simulated fields. Right: 
The residual between the non-Gaussian Euler characteristic and the Gaussian approximation. The solid line and the shaded band are 
the same as in the left panel, minus the Gaussian contribution. While the solid line is obtained by using the cumulants measured in the 
simulations, the dashed line uses the cumulants as predicted by the perturbation theory. 



with x ~ °" n - The first orders of this expansion in powers of a are (see Appendix E): 

= (2 1 (q 2 I 1 )+4(xI 2 })Ho(u)-(j 2 (xq 2 } + y(x 2 h})H 2 (u) + ^(x 3 )H 4 (u), (39) 
X (2 V) = (yv? 4 >c + 2 7v z<Z 2 J 1 ) c + 2{x z h} c ^H 1 (u)-(^(xV)c + ^{x z J x )^j H 3 {v) + ^{x%H 6 { V ) 

2 

- {^{xq 2 )(q\h) + {q 2 J 1 ){x 2 J 1 )+4{xq 2 ){xh}) H t {p) - Q{xq 2 )(x 3 ) + ^(x 3 )^ 2 ^)) H 5 {u) + ^ (x 3 ) 2 H 7 (v) 
+ (~f 2 (xq 2 ) 2 + ^(x 3 )(q 2 J 1 )+ 1 (xq 2 )(x 2 J 1 ) + \{x 2 J,) 2 + |<£ 3 )<:rJ 2 )) H 3 (u), (40) 

while the a 3 correction, X (")> i s given m Appendix E by equation (E15). Figure 1 compares the first-order non-Gaussian 
prediction to simulations of non-Gaussian density field developed through non-linear gravitational instability. As described in 
detail in section 4.2.3, we performed 25 scale-invariant (n s = —1) 256 N-body simulations which are then cut into 25 slices 
each to form 2D fields. The cubic cumulants required to form the non-Gaussian prediction (equation 39) are then measured 
from these slices and the 2D Euler characteristic is determined using the suite of algorithms described in Appendix G. The 
agreement between the measured Euler characteristic and the theory is remarkable. Perhaps the most striking fact is that the 
non-Gaussian correction indeed arises as a combination of low order Hermite polynomials. Alternatively, in this model the 
cumulants can be predicted by perturbation theory (PT, see section 4.2.1). The right panel in Figure 1 demonstrates that 
both approaches give consistent results 7 . 



3.2.3 Non-Gaussian total and differential extrema counts in 2D 

Extrema counts, especially that of the maxima of the field, also have long application to cosmology (e.g. Bardeen et al. 
1986), however theoretical development have been mostly restricted to Gaussian fields. Given the knowledge of the JPDF 
(equation (11), non-Gaussian extrema can be computed in the similar way as the Euler characteristic (Pogosyan, Pichon & 
Gay 2011). For example, the density number of all critical points above the threshold v is 

n(v) = n max (v) + n sad (u) + n min (v) = — s / dJi/ r AC, j dJ 2 P e xt(Cj Ji, J2)\h\ , (41) 




the difference from equation (35) only involve taking the absolute value of the Hessian determinant I 2 . Different types of critical 
points can be considered separately by restraining the integration domain in the 3\-Ji plane to ensure the appropriate signs 
for the eigenvalues. Specifically, for maxima the integration range is Ji G [-oo,0], J 2 G [0, J 2 ], for minima it is Ji G [0,oo], 
J2 G [0, J 2 ], and for saddle points Ji G [—00, 00], J 2 G [J 2 , 00]. 



7 For illustrative purposes we restrict the comparison of our all-order expansion (37) to simulations to first order non-Gaussian corrections 
as they can be predicted from first order PT. The second order terms have been studied in the perturbation theory in Matsubara (2010). 
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The total density number of extrema is obtained by integrating over threshold v. The effect of the first (cubic) non- 
Gaussian correction on extrema density in 2D is given by 



1 | f8(g 2 Ji)-5(J 1 3 ) + 6(JiJ 2 ) 



8V3-KRZ 
1 



+ 







1 f8(g 2 Ji)-5(J 1 3 > + 6(JiJ 2 ) 



(42) 
(43) 
(44) 



The total number of saddles, as well as the total of all the extremal points are preserved in the first order (the latter following 
for the former, as n max — n sa( j + ro m in is the Euler characteristic of the manifold and is thus constant). 

The differential distribution of the number of extrema with the threshold can be also found analytically in 2D: 



9n ma 
dv 



<9n min 



27ri?: 



exp 



1 + erf 



+ 



exp 



y/2n(l - 2/3 7 2 )i? 2 * V 6 - 4 7 
1 



dv 



dv 



2-kR; 



exp 



1 - erf 



yV 


V2(l-7 2 ) 


3v 2 \ 


1 


_4 7 2 J 




yv 





K\ (v, y) H — — exp , 

H ' ^{l-^)Rl P V 2(1 - 7 2 ) 



K 3 (v, 7) 



1 + erf 



yv 



v/2(l-7 2 ) 



V2(l- 7 2 )(3 - 2 7 2 ) 
1 



K 2 (y, 7), 



(45) 



Jfi(i/,7) 



V^(l - 7 2 )i?i 



exp 



2(1 - 7 2 ) 



#3(1', 7) 



H , — exp , 

V2tt(1 - 2/3 7 2 )i? 2 V 6-4 7 ' 



3^ 



^(1 - 2/3 7 2 )i?; 



exp 



6 - 4 7 5 



1 - erf 



K 2 (v,j), 



yv 



V2(l- 7 2 )(3-2 7 2) 



K 2 (v, y) , 



(46) 
(47) 



where K\, K 2 , K3 are polynomials in v with coefficients given by the cumulants and 7. Since 



9n max _^ 3n min 9n sad 



9i/ 



2-kRI 



exp — - I K x {v,y) 



(48) 



we immediately find that K\ is related to the series that appear in the expression for the Euler characteristic, equations (37) 
and (39). Specifically, K\ reads 



Ki(v,y) = 



8tt 



H 2 (v) + (H (gVr) + i (xJ, 2 ) - i <xj 2 >) - (<xg 2 ) + i (x 2 ^)) ff s M + J <* 3 > H 5 (v) 



(49) 



In the limit v 3> 1, the 3 terms are dominated by their respective exponentials. Since 2(1 - 7 2 ) and (6 - 4 7 2 )/3 are smaller 
than 1, the K 2 and K3 terms decrease more quickly than the K\ term. Thus in the limit of large v the expression of the 
density of maxima is dominated by the K\ term and coincides with the expression of the Euler characteristic, as expected. 

To write the remaining K 2 and K3 more compactly, we introduce two types of the scaled Hermite polynomials H J (v, k) = 
k n H n (v/k) and H^(v,k) = k~ n H„ (v/k), where the choice of k will be dictated by the exponential prefactor in the corre- 
spondent term. The polynomial K 2 is given by the straightforward expansion in Hn(v, ^/l — 2/37 2 ) 

1 



K 2 (v,y) = 



(W) + \l (q 2 -h) + I (xJi 2 ) - I (x,J 2 ) + | 7 ( Ji 3 ) - 2 -y (Ji J 2 >) (v, x/1 - 2/3 7 2 ) 
\ (V ) + 2y(x 2 .h) + ^y 2 (xJi 2 ) + \y 2 (xJ 2 ) + ^ 7 3 (Ji") + ^ (J1J2)) «F (v, y/l - 2/3 7 2 ) 



(50) 



The knowledge of K\ and K 2 is sufficient to describe the differential number density of the sum of extrema of all the types 



dn 



t)v 



2-kRI 



exp 



K±{v,y) + 



exp , 

^2n(l - 2/3 7 2 )i? 2 V 6-4 7 2 



3v 2 



K 2 (v,y). (51) 



The remaining term, Kz(v) is the most complicated. It is the simplest to express it as the expansion in W^(v, \fl — 7 2 ) but 
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Figure 2. Left: differential distribution of maxima in 2D slices of 3D density field evolved to <r = 0.1 for several values of the initial 
spectral index {red: n B = —1.8, green: n B = —1.5, blue: n B = —1, purple: n B = 0). Right: Measured maxima (blue), saddle point (green) 
and minima (red) distributions for n B = — 1 versus Gaussian (dashed lines) and the first order non-Gaussian (solid lines) predictions. 
The shaded region corresponds to the 3-sigma dispersion within the simulation set. 



the expression is still quite long 
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(3 
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7(3 
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2 7 2 ) 
3 M) 



(x*Jl) 



(52) 



Figure 2 (left panel) shows samples of this distribution for several values of the spectral parameter 7 according to perturbation 
theory (see Section 4.2.2). Figure 2 (right panel) shows the very good agreement between these predictions and the measure- 
ments of the corresponding number densities in simulations (see Section 4.2.3). In 2D, small non-Gaussianities has little effect 
on the saddle points, slightly shifting the differential distribution towards the lower field thresholds, while keeping the total 
number of the saddles constant to first order. In contrast, there are more significant effects on the maxima and minima. The 
maxima distribution becomes broader, developing high v excess relative to the Gaussian formula. The minima distribution 
behaves oppositely. It sharpens around its peak at v ~ —1.5, while at the same time losing the lowest minima. This reflects 
the intuitive picture of non-linear gravitational clustering developing high density peaks, while stretching and smoothing low 
v voids. 



& 2 Xi& 3 x k iP(v, it, x k i)S D (x 2 ) |A 2 | , (53) 



3.2.4 The differential length of the skeleton for non-Gaussian fields in 2D 

The skeleton is a tool used to extract the filamentary structure of a field (Novikov, Colombi & Dore 2003). The theory of the 
skeleton for Gaussian fields has been presented in Pogosyan et al. (2009), which demonstrated that its differential length can 
be approximately (in the so called stiff approximation) computed as 

dL skcl _ 1 
dv ~ R 

where the Xi are the components of the gradient in the Hessian eigenframe. To separate the part of the skeleton that maps 
closely the high-field value filamentary structures from the other critical lines, in Pogosyan et al. (2009) the notion of "primary" 
skeleton was introduced, defined by the constraint 

Ai + A 2 < 0. (54) 

The skeleton presents an insightful case study of how our invariant formalism can be extended to the cases in which the 
independent rotation symmetry of the gradient and the Hessian is broken. The skeleton definition couples the directions of 
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the gradient and the Hessian: it singles out the direction associated with the highest eigenvalue of the Hessian (Ai) which 
is the local direction of the skeleton. While the components of the gradient orthogonal to the skeleton's direction are zero, 
their parallel component should be integrated over. This does not break the isotropy for the Hessian-related part of the 
expansion since the integrand only contains A2 which is a rotation invariant of the Hessian. However it means that we cannot 
independently average over the directions of the gradient and the expansion has to treat its components separately. 

We will first consider that the statistical properties of the components of the gradient are the same in the Hessian 
eigenframe as in any frame. Appendix C discusses this assumption and shows how to avoid it and that it only has a secondary 
impact on the results. In practice the gradient part of the Gaussian PDF is 

G(xi, X2) = — expf— x\ — £2). (55) 

7T 

_ 2 

The corresponding orthogonal polynomials needed to build an expansion around the Gaussian kernel x are the 

"physicists'" Hermite polynomials H% hys (x) = 2 5 H n {y/2x). Their orthogonality relation is 

/ + o° 1 2 m 

-=e~ x 2? H n (V2x)2* H m (V2x)dx = n!5 n , m . (56) 
-00 V 71 " 

The Gram-Charlier expansion is then 

P{(,x 1 ,x 2 ,J 1 ,J 2 ) =Gll+J2 ^2 a a ! b \ l\ ]\ {^^^^^oc Hl (C) H a (V2x 1 )H b (V2x 2 )H k {J x )Li (J 2 ) 

\ n — 3 i,a,b,k,l 

with the normalization of the Gram-Charlier coefficients again such as to have the highest-degree term to coincide with the 
corresponding moment: 

(H> = -4(-ff»(v^xi)>. (57) 
2 2 

Using the variable x instead of £ allows to introduce a threshold more easily. By definition of £ = (x + yJi)/\/1 — J 2 , the 
Jacobian of the transformation is d£d Ji = 1/ yl — 7 2 did Ji and the Gaussian PDF becomes 

G(x,xi,x 2 , Ji, J2) = / = exp (- ^ J" 7Jl J -X! 2 -x 2 2 - ^ J2] ■ (58) 

2-Ky'l— Y \ — 7 ) 1 J 

The expression of the non-Gaussian PDF in these variables is then 

/ 00 i + a + t> + fc + 2Z — n 

P(x,X 1 ,X 2 ,J 1 ,J 2 ) = G(x,X 1 ,X 2 ,J 1 ,J 2 )\ 1 + E 

\ n— 3 i,a,b,k,l 

To express the integrand | A2 1 via the Hessian invariants and determine the appropriate limits of integration we note that the 
skeleton constraint, Ji = Ai + A2 < 0, implies that A2 < 0, since eigenvalues are sorted, Ai > A2 . The relation between the 
integrand and the invariants is thus 

\M = l(VJ2-Ji) , (60) 
and the differential length of the skeleton is given by 

dv = 2R J d< ^ 2 J J J dx 2 5 D (x 2 ) - JiJ P{x = v, xi, x 2 , J\, Ji). (61) 

The integration over x 2 is trivial while the integration over x\ corresponds to the projection on H§(\/2x\) that selects only 
the a — term: 

orskcl 1 r+00 rO . I 

00 i + b+k + 2 l = n . , / T \ \ 

E E ^(^J^^f^jw.WX.w). (62) 



The result of the remaining integrations takes the form 
dL skcl 1 



_ -e 2 

9^ N/27T-R, 



Aexp < _iv \ +B / 1 + erf f - 



2(l-7 2 )/ V V\ / 2(l-7 2 ) 
In the Gaussian case, the coefficients A et _B were given in Pogosyan et al. (2009): ^ (0) = sj\ - ~i 2 / (2y/2n) and B (0) 



(63) 
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Figure 3. Left; differential distribution of skeleton length in 2D slices of 3D density field evolved to a = 0.1 for several values of the 
initial spectral index {red: n B = —1.8, green: n B = —1.5, blue: n B = —1, purple: n B = 0). The dotted line is the Gaussian approximation, 
the solid line is the first order non-Gaussian prediction. Right: Measured length distribution of the skeleton (blue) and the anti-skeleton 
(red) for n B = —1 versus Gaussian (dashed lines) and the first order non-Gaussian (solid lines) predictions. The shaded region corresponds 
to the 3-sigma dispersion within the simulation set. 



l/(8i/n) (-y^TT + 2-yv) . At the a order, the non-Gaussian correction are 

a.. _ \ ' - " "0F(-27 3 + 7 5 ) + 2(37 2 - 57 4 + 27 6 )i/ 1 ( I ,)-0F(37 - 37 3 +7 5 )^2M + 2(l-7 2 ) 3 // 3 M^3 



A 



2\/27r 



12(1 - 7 2 ) 3 



+ 



+2(2^-3^ + 1 5 )H 1 (u)-^H 2 (v), s 



4(1 -7 2 ) 3 



x Jl) + 



-V^7 3 - 1Hi(y) - V^lH 2 (v) 



4(1 - 7 2 ) 3 



- 2 7 2 ) - 27^(1/) - V^j 2 H 2 (v) 3 , ^7 / t \ 
H —i (Ji ) - TR T2\( xj2 l 



12(1 - 7 2 ) 3 



-(J1J2) 



4(1 ~ 7 2 ) 4(1 -7 2 ) 



(64) 
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(65) 



Figure 3 (left panel) displays the variation of this differential length w.r.t. the underlying power-index, n s . The dynamical 
evolution tends to lower the total length of the skeleton. This corresponds to filaments merging. At the same time, at the 
very high density tail v > 2 notice an increase in the length of the skeleton per unit volume, which reflects the compression 
of the high-density regions. These predictions can easily be adapted to the anti-skeleton that maps the troughs of the field 
by simply replacing both fields and cumulants with their opposite, i.e. replacing v with —v and the odd moments by their 
opposite. Figure 3 (right panel) shows a reasonable agreement between these predictions and the measured differential length. 
The discrepancies between the predictions and the measurements could come from the approximation used to predict the 
differential length of the skeleton (the "stiff approximation", as explained in Pogosyan et al. 2009) or from differences between 
the analytical definition of the skeleton and its numerical implementation (Sousbie, Colombi & Pichon 2009). Note that it 
is the total skeleton length that is most sensitive to a possible mismatch of the exact definition of the "primary" skeleton in 
the theory versus the techniques used in numerical simulations, while the differential behaviour at high threshold values is 
much more robust, as was pointed out in Pogosyan et al. (2009). Accordingly, when comparing equation (63) to the data in 
the right panel of Figure 3 we correct the overall amplitude by a factor Ltot, measured/^tot, predicted = 0.23 for n B = —1. As 
non-Gaussianities distort the skeleton statistics without much changing its total length, this empirical correction does not 
prevent us from estimating a from the shape of the differential length, as is seen in Figure 3. Incidently, in 3D, we found the 
amplitude to be correctly predicted (modulo some change in the condition for the definition of ridges), so no such correction 
is required (see below Sec. 3.3.4). 
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3.3 Non-Gaussian critical sets in three dimensions 

In 3D dimensions the basis of our consideration is equation (24). The expressions for the critical sets are somewhat more 
complex, while cosmology provides a context in which we can directly predict the value of the relevant cumulants within the 
framework of gravitational perturbation theory (see Sec. 4.2.1). 



n— 3 i,j 

For instance, to first order in the non-Gaussian expansion, we have 



(66) 



3.3.1 The area of 3D isosurfaces 

Let us again start with a 3D Minkowski functional, the surface, A(v), of the x = v isocontour, which is given by 
A(u) = — dx Aq 2 P(x,q 2 )8°(x-v)q = — \ Aq 2 P{v,q 2 )q 

n J-oo Jo U Jo 

where the joint PDF of the field and its gradient, P(u, q 2 ) follows from equation (24) after integration over the J^'s: 

iw> - ^v?«p (4 - 1/) (i + 1 T" < iV 'L ».w*r <§«■>) ^ w 

\ ^ — 3 i ,j / 

Performing integration using J °° exp(-3/2g 2 )g 2 L^ 1/2) (\q 2 )dq 2 = -2r(j - \)/(V^T(j + 1)) yields 



Appendix E2 demonstrates how to re-order expansion (68) as an Edgeworth expansion. 



3.3.2 The Euler characteristic of non-Gaussian 3D fields 

To compute the Euler characteristic in 3D, following the 2D case, one averages — 13 = — (Ji 3 -3JiJ 2 + 2J 3 )/27 over the full 
range of the second derivatives: 



xiy) 



,3/2 

1 r 00 r 00 r 00 f 2 

•jjgj dJi J^jdC J dJ 2 y 3/ dJ 3 P oxt (C, Ji, Ja, Js)/ 3 



(70) 



As shown in Appendix A2, the resulting expression is 



xH = ^Erfc( J^) x (_oo) + 
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00 i + 2j + fc + 3=n . 1 vfe+l 

+ E E 
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I]lF \ C ? Jl J3 /oc (1 - 7 ) 7 M 



(71) 



As in 2D, in the rare event limit, v 2> 1 or v <C 1, equation (71) also represents resp. the general differential number of 
maxima or minima. When equation (71) is reordered in a, the first corrections are (see Appendix E) : 

X (1) = + 97^/2)) Hi(u) - (§7 3 <*<? 2 > + -tV*)) H 3 {u) + p 3 (x 3 )H B {u) , (72) 
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V V 

Figure 4. Left; Euler characteristic in 3D for a mildly non-Gaussian field (cr = 0.1). Dashed line: Gaussian prediction, solid line: first 
order prediction. The shaded band corresponds to the 3-sigma dispersion for the measurement of the Euler characteristic in the simulated 
fields. Right: The residual between the non-Gaussian Euler characteristic and the Gaussian approximation. The solid line and the shaded 
band are the same as in the left panel, minus the Gaussian contribution. 



and 

X (2) = " (f l(q 2 h) c + 27<zl 3 >c) H {v) + 07V>c + f^W Ji)c + \l(xh)^ H 2 {v) 

- (j7W)c + ^7VJi>c) Hi{y) + ~y 3 (x%H 6 ( V ) + (^7<<M> 2 + y (q 2 Jt)(xI 2 ) + y <^ 2 )</ 3 )) H (u) 
- (f l 2 {xq 2 ){q 2 Ji) + \^{q 2 J.){x 2 Ji) + 2 ^{xq 2 ){xh) + l^J^ixh) + ^ 3 / 3 )) H 2 (v) 
+ (f 7W> 2 + -^{x^iq 2 ^) + lj 2 (xq 2 ){x 2 J,) + -^{x 2 ^} 2 + ^(x 8 )^)) H t {u) 

- Q 7 3 <sg 2 >(x 3 ) + \~t 2 {x z ){x 2 J^ H 6 (u) + ^{x^H^v). (73) 

Figure 4 shows that the first order expansion provides a good match to the measured 3D Euler characteristic, matching 
closely the functional form of the non-Gaussian correction, antisymmetric in 3D. Measurements are performed by computing 
the alternating sum of the counts of critical points. 

In cosmological applications we do not expect our Universe to have complex topology at the large-scale structure scales, 
if at all, and, thus \{— oo) = 0. This matches exactly (J3) = from equation (70), i.e, the following relation between the 
moments of the 3D Hessian invariants holds 

(J 1 3 )-3(JiJ 2 ) + 2(J 3 ) =0 . (74) 



3.3.3 Non-Gaussian total and differential extrerna counts in 3D 
The total number of critical points in 3D is given by 



1 r+00 r+00 r+00 rJ2 3 ^ 2 

next = ^3/ dJi / dC / dJ 2 / dJ 3 P cxt «, Ji, J 2 , J 3 )\h\. (75) 

Kt J„ oc J " +tJi J y_j 2 3/2 



For a specific set of extrerna, sign constraints must be imposed on the sorted eigenvalues Ai > A 2 > A3 which establish the 

non-overlapping ranges of integration in (Ji, J 2 , J 3 ) space, defining the maxima ( ), the minima (+ + +) and the saddle 

points of "filamentary" (H ) and "pancake" (+H — ) types. While J 2 continues to span unrestricted range of positive values, 

for each fixed J 2 value, the boundaries between the different extrerna types are set in the (Ji, J 3 ) subspace by the cubic curve 

o 3/2 3/2 

Is = 0, i.e 2J 3 = — Jf + 3Ji J 2 , and by the limits — J 2 < J 3 < J 2 . In Figure 5, the ranges of integration are shown for each 
extrerna type. 

The total number of extrerna can be found analytically at all orders of the non-Gaussian expansion. Here we give only 



© 0000 RAS, MNRAS 000, 000-000 



16 C. Gay, C. Pichon and D. Pogosyan 




Figure 6. Left: Sample differential distributions of maxima (right group of curves) and filament-type saddle points (left group) in 
3D for several values of the spectral index {red: n s = —2.5, green: n B = — 2, blue: n B = —1, purple: n B = 0). Dotted lines are the 
Gaussian approximation, solid lines are the non-Gaussian prediction for a = 0.1 computed via Monte-Carlo integration. Right: Measured 
differential distribution of critical points in 3D for a mildly non-Gaussian field from n s = —1 simulations at a = 0.1 versus the first order 
non-Gaussian prediction (solid) and the Gaussian approximation (dashed). 



the first corrections that read 

2<VI5t18VI0 5^/5 /. 2 . 8 , 3 . 10 

29^ T 18^ 5^5 (/2 7 \ 8 /r3\ , 10 /r r A 

[} q Jl '~ 21 ' 1 ) + ^t^ JiJ2 ) ) 



(76) 



18007r 2 Ji2 24n 2 V6^Rl 

signifying that the first non-Gaussian corrections are equal in magnitude for all extrema types acting in the same direction 
for the maxima and the "filamentary" saddle points, and opposite to that for the "pancake" saddles and the minima. 

For the differential distributions of critical points, no analytical expression is known even for Gaussian fields and we thus 
proceed numerically. We use a semi Monte-Carlo integration: we draw values of the field and its derivatives according to the 
corresponding Gaussian distribution and then integrate the integrand times the Edgeworth correction over the domain of 
integration consistent with the required sign for the eigenvalues 8 . The left panel of Figure 6 shows the results of the numerical 
integration for non-Gaussian maxima and filamentary saddles of the density fields that are obtained by non-linear gravitational 
evolution from several scale-invariant initial conditions with power-law spectral index n s £ [—2.5,0] to a — 0.1. The numerical 
results for non Gaussian distribution of maxima are accurately fitted by the expressions given in Appendix B. 

In the right panel of Figure 6 the first order perturbation correction is compared to the actual measurements of the 



This technique is inspired by importance sampling (Robert & Casella 1998). 
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differential number densities of extremal points for n s = — 1 scale invariant power spectrum. The theoretical predictions fit 
the measurements very well for a — 0.1, correctly capturing the main qualitative features - the development of an excess count 
for the high amplitude peaks and filamentary saddles, with a simultaneous decrease in the number of peaks and filamentary 
saddles for average contrasts (y ~ 2 for peaks and v ~ 1 for saddles of this type) and the opposite trend for minima and 
"pancake" saddles, namely, depopulation of the underdense tails of their distribution compensated by an increase in numbers 
at contrasts near the mean (y ~ —2 for minima and v ~ — 1 for "pancake" saddles). We note that according to equation (76) 
and the values of the cumulants from Table Fl, the total number of extrema has changed very insignificantly at a — 0.1, 
by just over 2% for maxima and minima and ~ 0.07% for the saddle points. This shows that the differential distribution of 
extrema counts is a much more sensitive tool to probe non-Gaussianities. 



3.3.4 The differential length of the skeleton for non-Gaussian fields in 3D 

In 3D, the length of the skeleton is given in the stiff approximation by: 

aL skci ^ l 

dv ~ M 

with the constraint Ai + A2 + A3 < 0. As in 2D, in order to compute the differential length of the skeleton we need to separate 
the component of the gradient aligned with the highest eigenvalue of the Hessian from the others. Let Q be the normalized 
norm of the projection of the gradient in the hyperplane orthogonal to this direction: 



& Xi& x k iP{v,Xi,x kl )8 (x 2 )8 (x 3 ) |A 2 A 3 | 



(77) 



" {q 2 - x 2 ) ■ 

In the Gaussian case, we have Q 2 = § (x 2 2 + xfj. The Gaussian PDF for q 2 is G(q 2 ) = 3 
for xi and Q 2 : 

' 3 



(78) 



„— qe 2« which leads to the PDF 



G(xi,Q 2 )dQ 2 = 



-Q' 



dQ" 



(79) 



After integration over X\, the Gram-Charlier expansion of the PDF is thus 

00 i+2j + k+Z — n 



P(x,Q ,Ji,J 2 ,J3) = G(x,Q , Ji,J%, J 3 ) 



i+E E 
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The differential length of the skeleton is given by 
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(80) 



(81) 



The integrand | A2 A3 1 is a non-trivial combination of the invariant variables: one would have to find the sorted roots of the 
third-degree characteristic polynomial A 3 — /1A 2 + I2X — I3 = 0. Furthermore Pogosyan et al. (2009) has shown that the 
differential length of the skeleton does not have an analytical expression in 3D. It is thus easier to use a numerical integration, 
as was described for the extrema. Note that we use Ai + A2 + A 3 < instead of Ai + A2 < (that was adopted in Pogosyan 
et al. (2009)) as the constraint for defining the skeleton; this choice yields a perfect match to the total length of the skeleton 
as measured in simulations (in Figure 7 the curves are not normalized by the total length). 

For n s £ [—2.5,0], the Gaussian and first order prediction for the differential skeleton length are given in Appendix B. 
Figure 7 shows the variation of a non Gaussian 3D skeleton with the spectral index, and illustrates the match between the 
predicted and the measured skeleton. 



4 APPLICATIONS: CMB PRIMORDIAL NON-GAUSSIANITIES, DARK ENERGY IN THE LSS 

With the full knowledge of the departure from Gaussianity of critical sets (resp. in 2D equations (37) for the Euler character- 
istic, (47) for extrema, (63) for the skeleton, in 3D, equations (71), (76) and (81)), let us know illustrate how these estimates 
can be used to address problems arising in cosmology. The first example uses 2D extrema counts to estimate primordial 
non-Gaussianities in synthetic CMB maps, while the second example estimates the density dispersion, <r(z) of dark matter 
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Figure 7. Left: differential length of the skeleton in 3D for several values of the spectral index {red: n B = —2.5, green: n B = — 2, blue 
n s = —1, purple: n s = 0). The dotted lines are the Gaussian case, the solid lines are the non-Gaussian prediction for a = 0.1. Right 
differential length of the skeleton (blue) in 3D for a mildly non-Gaussian field (<r = 0.1); solid line: Gaussian prediction; dashed line 
first order prediction. The match between the two curves is significantly improved, compared to the 2D skeleton. 



simulations from differential extrema counts, the Euler characteristic and skeleton length as a function of redshift. Note that 
alternative illustrations can be found in Pogosyan, Gay & Pichon (2009); Pogosyan, Pichon & Gay (2011), including models 
involving second order corrections. 



4.1 Primordial non-Gaussianities in the cosmic microwave background 

Let us generate sets of parameterized non-Gaussian maps using the package sky-ng-sim (Rocha et al. 2005) of HEALPix. In 
this so called harmonic model, the PDF, fr(T) of the pixels' temperature, T is given by 



/ T (T) = exp 



2al 



i T 

cgcjH., 

i=0 



(82) 



where the c; are normalization constants. In this paper, we use nside=2048, ^ max = 4096, n = 2, cr = 1, q = and vary 
oi\ and «2 on a 50 x 50 polar grid in the range [0.1,0.8] subject to a\ + a 2 < 1. For each set of maps, we compute its 
derivatives, and arithmetically average the corresponding cumulants, using the code map2cum (Pogosyan, Pichon & Gay 2011) 
for which the invariant variables Ji and J 2 on a sphere are defined via the mixed tensor of covariant derivatives Ji = x-^' % 
and J 2 = J\ — 4 \ j . The differential counts are then evaluated for a range of threshold, v G [—5, 5]. For our reference map, 
a.\ = a 2 = 0.6, the number of extrema is computed by the procedure hotspot (Gorski et al. 2005). Figure 8 illustrates our 
hability to recover input model from a given set of maxima and minima extracted from a given realization of the sky. Here the 
residual departure from the Gaussian number count is shown. Indeed the maximum likelihood solution, solid lines, corresponds 
to a very good match to the measured residual differential number counts. Note that this class of model has an interesting 
feature that (s 3 ) is significantly smaller than say, (J 3 ). Indeed for this model, (a; 3 ) = — 0.076, (x 2 Ji) = 0.072, (q 2 x) = 
-0.091, (xJf) = -0.15, (xJ 2 ) = -0.087, (q 2 Ji) = 0.054, (j?) = 0.31, (Ji J 2 ) = 0.023. This simple minded experiment should 
be improved to account for, e.g. more realistic and physically motivated non-Gaussian maps (such as /nl models), realistic 
galactic cuts and hit maps, variable instrumental response. The main challenge is to model incompleteness in observed samples 
in order to weight accordingly the range of relevant thresholds. 



4.2 Constraining the dark energy equation of state via the geometry of the LSS 

At supercluster scales where non-linearity is mild, we can expect the non-Gaussianity of the matter density to be also mild, 
but still essential for a quantitative understanding of the filamentary Cosmic Web in-between the galaxy clusters. Indeed, 
e.g. Figure 6 shows that the interplay of gravity and (possibly accelerated) expansion will induce some level of distortion in 
the PDF of the differential counts. Importantly, this distortion is not sensitive to the overall number-count as the PDF are 
normalized. Fitting measured counts to the distorted theoretical expectation should therefore allow us to robustly estimate 
some features of the underlying dark energy equation of state. 
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Likelihood for minima 6n Likelihood for maxima 




Figure 8. Middle panel: the predicted (solid line) excess number of maxima (bottom), and minima (top) in Au = 0.25 bins as a function 
of the threshold, u, on top of the measured count (dashed line) from a single realization full-sky nside=1024 HEALPix map (histogram) . 
The temperature field is smoothed with the Gaussian filter of fO arcmin FWHM, resulting in i?* pb f 1 arcmin RJ 3 pixels. Left and right 
panels: likelihood contours in the parameter space of the Harmonic Oscillator model of non-Gaussianity given by equation (82); the input 
model corresponds to cti = 0.6, «2 = 0.6, represented by a green dot in the likelihood map. The five contours correspond to 1,2,3,4 
and 5 sigma likelihood contours, while the blue and red dot mark the most likely estimate for maxima and minima respectively. Hence 
different models of non-Gaussianity can be well distinguished by their effects on differential extrema counts. 



4-2.1 Extension of the gravitational Perturbation Theory to critical set cumulants 

When the field is the 3D density field of the large-scale structure of the Universe and for Gaussian initial conditions, the non- 
Gaussianities appear because of the non-linearities of gravitational clustering. As mentioned in the introduction, equations (71), 
(76) and (81) have in common that they can be re-arranged as 

£ R {u,z) = Y,Cn{u,R)o^(R)D(z) n , (83) 

n 

where a 2 is today's variance of the cosmic density field at the corresponding smoothing scale, R, and D(z) is the growth 
factor as a function of redshift. For instance, for the Euler characteristic, CioqD(z) is given by equation (72) while C2C 2 D(z) 2 
is given by equation (73). When gravity is generating mildly non-Gaussianities, we will show that it is possible to predict 
the emerging statistical properties of the field using a perturbative theory. As we will demonstrate below, all the cumulants 
involved in the Edgeworth expansion entering the C„s may then be predicted from first principle. While the geometrical 
analysis of catalogs, say using codes like DisPerSE (Sousbie, Pichon & Kawahara 2011) allow us to measure Br(v,z), and 
provided that i) we consider large enough scales so that equation (83) converges, and ii) that gravity is the driving cause 
of non-Gaussianity, equation (83) can be inverted for D(z). In turn, measuring D(z) below z ~ 1 probes the dark energy 
equation of state. One of the difficulties in implementing our statistics to observed surveys is that a is not known, hence we 
cannot normalize v = p/o~. Following Gott (1988), a work-around is to introduce i/f, the density contrast corresponding to a 
given filling factor, /, instead of v which are related via the identity 

/■OO /"OO 

/ = / duP(u) = / duG(u) , (84) 

where f(v) is given by equation (29) and G stands for the corresponding Gaussian PDF, so that £n(u[uf], z) is a measurable 
quantity encoding the departure from Gaussianity which is not already encoded in P{v). In practice we will not make use of 
equation (84) below and we will assume perfect knowledge of u in these experiments. 



4-2.2 Prediction of the 3rd order cumulants in Perturbation Theory 

In the formalism of perturbation theory (Bernardeau et al. 2002), it is easy to see that the p-th cumulant of the field is 
of order a p ~ 2 , if the field is normalized to a constant variance. Using this prescription, one can re-sum the Gram-Charlier 
expansion to the so-called Edgeworth expansion. The Edgeworth expansion is an expansion in a and its order is thus physically 
motivated. In practice it means its convergence properties can be justified on physical grounds, that it can thus be truncated 
more easily (Blinnikov & Moessner 1998). Appendix D shows that the a term from the Edgeworth expansion only comes 
from the third order term of the Gram-Charlier expansion, while the a 2 term comes from both the fourth and part of the 
sixth order terms. We will therefore restrict our predictions of the involved cumulants to third order. At this order, moments, 
cumulants and Gram-Charlier coefficient coincide. For the next order in a, one would have to extend the predictions to quartic 
cumulants and to decompose sixth order Gram-Charlier coefficients in cumulants (see Appendix E). For third order terms, the 
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calculation of all the cumulants is similar to the well-known prediction of the skewness (Peebles 1980; Fry 1984; Bernardeau 
1993; Bernardeau et al. 2002). We will only consider an Einstein-de Sitter (EdS) universe since other cosmologies can be 
described as corrections to this prediction. 

For a three dimensional field, the PT cumulants involving the field and its derivatives can be written as: 

where d™ is the n-th derivative with respect to the i-th coordinate, and 8r is the field smoothed on scale R. Using gravitational 
perturbation theory, we show in Appendix Fl that Cr obeys 

C R = 2 j d 3 k!d 3 k 2 tFo,/s,7(ki,ka) + Jft 7 ,a(ki,ka) + T-,,*,pQ*i t *a)] P(k 1 )P(k 2 )W(k 1 R)W(k 2 R)W(\k 1 + k 2 \R) , (85) 

where the so-called generalized geometric factors (Bernardeau et al. 2002), -7^,^,7 (ki, k 2 ) (here a is the short hand for the 
triplet (01,02,03), etc. ) read: 

-F a ,^ 7 (k!,k 2 ) = F2(k 1 ,k 2 )e(k 1 ,k 2 )i Q1+a2+Q3 i' 3l+ ' 32+ ' 33 (_ip+T2+73 j with F 2 (k!,k 2 ) = | + + ^ (kl 2 k2 / , 

e(kl ,k 2 ) = ( k W) Q1 (kf^ H 31 )" 3 K ] )' 81 K)' 82 K 1 )" 3 K ] +4 11 ) 71 ((kf 1 + k^) 72 (kf + k^) 73 . (86) 

The smoothing of the field over the scale R is taken care of via the window function, W(kR). As explained in Appendix Fl, 
for a Gaussian filtering and a scale invariant power spectrum, these integrals can be completed and the resulting cumulants 
are given explicitly in terms of the smoothing length and the power index in Appendix F2. For instance, for a n E — —3 power 
law index, we generalize the classical result 

/ 3\ 34 . 2 v 34 2 34 2 , 2 . 34 2 , 2 \ 34 2 

{x) = —, to (aaBi 2 ) = y^ T gjj, (xxuXm) = y j-^ , (asm) = y and ( XXl ) = ——. (87) 

In Section 3.2 we also use the perturbation theory to predict 2D statistics. Since our 2D fields are the slices of the 
gravitationally evolved 3D density field, the cumulants in field variables that involve only the derivatives restricted to the 
slice are the same as predicted by the 3D perturbations theory and are given in Table Fl. However, to apply to 2D slices, one 
needs to combine these cumulants of the field variables into 2D invariants. 



4-2.3 Comparison to measurements of the critical set cumulants in Monte Carlo simulations 

Measurements are carried on sets of simulations produced using Gadget-2 (Springel 2005) code and scale-invariant initial 
conditions generated with Mpgrafic (Prunet et al. 2008). We use 256 3 particles and 10, 25 and 10 realizations for respectively 
n s = 0, -1 and -2. To be consistent with the perturbation theory used, we set an Einstein-de Sitter cosmology (fl m = 1, Qa = 0). 
The 2D results are obtained by cutting these 3D simulations into 25 slices which are averaged along the shortest dimension 
to give 2D fields. Following Colombi, Pogosyan & Souradeep (2000), we benefit from the scale-invariance of our simulations 
to combine them. For each simulation, we compute the density field with a cloud-in-cell procedure and we smooth on lengths 
of R = 5, 10, 15, 20 and 25 pixels. To eliminate the snapshots contaminated by initial transients or incorrect large-scale 
mode couplings, we compute the correlation length lo, defined by o 2 (lo) = 1. We only keep the snapshots verifying (Colombi, 
Bouchet & Hernquist 1996): 

] v T 7 3«'o«i, 

where L is the size of the box and iV the number of particles. The first condition ensures that the grid effects have relaxed, 
while the second condition limits the coupling of the modes modified by the finite size of the box. In practice we choose 

„ L , 1 r 
< lo < -L. 



N i/3 ^ 20 

Even within a selected snapshot, all the scales cannot be trusted. The smoothing length has to be large to probe scales 
for which the discreteness of the simulation is not important, but small enough compared to the size of the box. A detailed 
study shows that the smallest scale containing information is not given by the initial inter-particle distance, as one could 
naively think. Indeed the resolution is improved in clustered regions, where the density is higher than the density of the 
initial conditions. The pertinent characteristic length can be computed using the quantity N C (R) = n£(R), where 7? is the 
smoothing length, n the average number of particles per cell and £ the correlation function. N c characterizes the average 
number of neighbors (Balian & Schaeffer 1989). The associated characteristic length is then l c such as N c (l c ) = 1. We keep 
the smoothing lengths for which (Colombi, Bouchet & Hernquist 1996): 

lc < R < L. 

In practice, we choose 2l c < R < L/20. The rotational invariants of the field are computed for each field (using Fourier 
transforms to estimate the derivatives), and are averaged over realizations once the relevant cumulants are generated. The 
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(xJ! 2 )/* 


3.871 


3.75 ±0.06 


3.877 


4.1 ± 0.2 


3.847 


5.3 ±0.6 


{xJ2)/(T 


1.545 


1.54 ±0.02 


2.037 


2.06 ±0.05 


2.625 


3.1 ± 0.3 


{q*Jl)/v 


-1.335 


-1.28 ±0.02 


-1.157 


-1.28 ±0.08 


-0.941 


-1.5 ±0.2 




-4.644 


-4.50 ±0.08 


-4.412 


-4.8 ±0.3 


-3.960 


-7±1 


<JlJ 2 >/<T 


-0.679 


-0.65 ± 0.01 


-0.955 


-0.92 ±0.02 


-1.142 


-1.5 ±0.1 




1.304 


1.28 ±0.03 


0.774 


1.0 ±0.1 


0.267 


1.0 ±0.3 



Table 2. Predicted and measured 3D third order invariant cumulants for three values of the power spectrum index. The method to 
estimate these cumulants in simulations is described in Sec. 4.2.3. The predictions follow from the coordinate cumulants computed in 
Table Fl given in Appendix F2. According to perturbation theory, the dimensionless third order cumulants are proportional to a. 





n B = 


n s = —1 


n s = -2 




n s = 


n s = —1 


n c = -2 




39 ±2 


53 ±8 


140 ± 60 


(x 2 J 2 )c/<T 2 


5.1 ± 0.3 


9.3 ±0.9 


22 ±6 


{Jl 2 J2)c/a 2 


3.6 ±0.2 


7.0 ±0.6 


21 ± 7 


{x^J^c/a 2 


-20 ± 1 


-26 ±3 


-45 ± 14 


(J2 2 )c/* 2 


4.4 ± 0.2 


7.4 ±0.8 


21 ± 8 


(X^c/a 2 


17 ±1 (17.5) 


23 ± 3 (21.9) 


40 ± 12 (30.4) 


(JlJ S )c/<T 2 


-1.6 ±0.1 


-2.0± 0.3 


-4± 2 


(Q 2 JS)c/v 2 


16 ±5 


17 ± 5 


16 ± 5 


(xJ^) c /a 2 


-30 ±2 


-40 ± 5 


-88 ± 36 


(q 2 J2)c/« 2 


13 ±5 


14 ± 5 


13 ±5 


(xJ 1 J 2 )c/cr 2 


-3.6 ±0.2 


-6.5 ±0.5 


-17 ± 5 


{xq^J^c/a 2 


-7.3 ± 0.4 


-9± 1 


-19 ± 7 


(xJ 3 ) c /* 2 


2.4 ± 0.2 


4.0 ±0.8 


7±4 


{x 2 q 2 ) c /a 2 


8.5 ±0.5 


12 ± 2 


24 ± 8 


{x 2 J{ 2 )c/* 2 


24 ± 1 


32 ±4 


62 ± 22 


(Q 4 )c/v 2 


6.8 ±0.4 


10 ± 1 


24 ± 9 



Table 3. Measured 3D 4th order cumulants as in Table 2. The corresponding predictions, which would involve extending second order 
PT to invariants is left to future work. Those for (x 4 ) c (from Lokas et al. 1994) are indicated in brackets. According to perturbation 
theory, the dimensionless third order cumulants are proportional to a 2 . 

results for the 3rd order moments of the field are given in Table Fl and are compared to predictions from first order perturbation 
theory. Table 2 re-expresses them in terms of invariants and Table 3 gives the 4th order cumulants. 

4-2.4 Growth factor measurement experiment: towards a possible dark energy probe? 

To demonstrate the potential of our non-Gaussian critical set statistical estimates to map the evolution of the growth factor 
D(z) (equations (71), (76) and (81)), let us discuss a very idealized topological dark energy probe experiment. For this purpose, 
we generate a cosmic sequence of a set of 25 n s = — 1 scale-invariant simulations with 256 3 particle smoothed over R — 10 
and 20 pixels and measure for each redshift the corresponding differential number of extrema, £^ xt (^, z) (resp. max, min, 
filament-saddle and pancake saddle), Euler characteristic £^ uler (^, z), and skeleton length, £ s ^ cl {v, z). We denote the set of 
such measurements as {£^(1/, z)}k, where k stands for the type of statistics and R the smoothing scale (c.f. equation 83). We 
then match these histograms to the predictions of the first order Gram Charlier PT model and deduce the best fit a for a 
range of redshifts, using for now a least-squares fit with uniform weights. This estimated "geometric" a is then plotted as a 
function of underlying a for the six critical set statistics (4 types of extrema, Euler characteristic and skeleton length 9 ) and 
shown in Figure 9. On average, a is a good match to a (when a 0.2), with a dispersion cr| which scales with the smoothing 
length R and the size of the survey L as 

4 = a k c [{L c /R c ){R/L)f 2 , (88) 

where the values cr* that were measured in our numerical survey (L c = 256 pixels) are given in Table 4 for two values of the 
smoothing scale, R c = 20 and R c = 10 pixels. The deviation at larger a seen on Figure 9 could be explained as an effect of 
the higher order non-Gaussianities and should be reduced by going beyond a first order description. This is corroborated by 
the observations that the critical sets which are least sensitive to the higher and, presumably, more non-linear thresholds (e.g. 

9 The skeleton length is only fitted on the range v < or v > 2 where the agreement between predictions and measurements is the best. 



R c = 20 


<7™ ax ~ 0.030 , 


<r™ in ~ 0.034 , 


CT sa,d-iU ^ Q Q49 ^ 


sad — pan nni; . 

<r c ~ 0.064 , 


o-= kcl ~ 0.033 , 


CT Eulor ^ o ^ 


R c = 10 


CT max ^ 0.012 , 


cr™ in ~ 0.010 , 


CT sa,d-iU ^ _ 017) 


sai l — pan n n-i ^ 

<r c ~ 0.017 , 




CT Eulor ^ 0.012 . 



Table 4. Measured dispersion around the mean value extracted from the scale invariant experiment described in Figure 9 (L c = 256 
pixels) for smoothing lengths R c = 20 and 10 pixels. 
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(a) Minima 




(b) Pancake-type saddle points 



(c) Filament-type saddle points 
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(d) Maxima 



(e) Euler characteristic 



(f) Skeleton 



Figure 9. Fitting <r using the Gram-Charlier and the perturbation theory results, for the minima distribution (top left), the pancake-type 
saddle points distribution (top middle), the filament-type saddle points distribution (top right), the maxima distribution (bottom left), 
the Euler characteristic (bottom middle) and the skeleton length (bottom right), to a set of 256 3 , n s = —1 scale-invariant simulation 
smoothed over R = 20 pixels (which corresponds to approximately 230 maxima) for the larger and lighter contours and R = 10 pixels 
(1800 maxima) for the small and darker ones. The skeleton is smoothed only for R = 20 to ensure sufficient smoothness. The shaded area 
corresponds to one standard deviation (computed from 25 simulations). Note that the critical sets which are most sensitive to the higher 
thresholds tend to depart from the non-Gaussian linear correction earlier, which is expected in the context of gravitational clustering. 



minima and saddle points) tend to depart from the non-Gaussian linear correction later. On the other hand, our statistics 
are less sensitive at very low cr's due to the universal behaviour of the geometrical descriptors in the Gaussian limit. Hence 
the optimal level for the recovery of the variance a is around a = 0.1, according to Figure 9. Thus, varying the smoothing 
scale with redshift, R(z), in such a way as to ensure a level of non-linearity near cr(z) as 0.1 may help the reconstruction of 
D(z). Note finally and importantly that the overall spread in the measured a reflects cosmic variance and could be reduced in 
proportion to the square root of the volume of the survey, following equation (88). For instance, a naive extrapolation of the 
cluster CFHTLS Wide counts 10 (Durret et al. 2011) to a full sky experiment with 40 bins of redshift below z — 1.15 would 
yield a noise-to-signal ratio of 1% per bin for D(z) using the Euler characteristic alone. 

Let us now briefly discuss how these very idealized measurements can be extrapolated to anticipate the accuracy of a 
dark energy experiment based on critical sets. A dark energy probe may directly attempt to estimate the so-called equation of 
state parameters (w a ,u>o) from these critical sets while relying on the cosmic model for the growth rate (Blake & Glazebrook 
2003), 



D(z\w ,w a ) = — - — H(a) 



da 



iH 3 (a') 



H 2 (a) = Hi 



— + i 2a exp 

n-s 



l+w(z') 
1 + z' 



dz' 



(89) 



with Q TO) J1a and Hq resp. the dark matter and dark energy densities and the Hubble constant at redshift z = 0, a = 1/(1 + 2) 
the expansion factor, and with the equation of state w(z) = wq + w a z/(l + z). Given equations (83) and (89), optimizing the 
probability of observing all counts at all redshifts with respect to (w a , Wo) yields a maximum likelihood estimate for the dark 



http:/ /www. astromatic.net/iip/wl. html 
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energy equation of state parameters 11 : 

(w ,w a ) = arg max < TT Poisson ( £«K, D(z\w a , w )], E* 4 ) > , (90) 

^Oi^a I V / I 

I 2,1/^ ,fc 1 

where the predicted "number-counts", £ %[vi, D(z\ w a , i£>o)] are given by equation (83), given equation (89) 12 . Here E^ ti is the 
measured k statistics in the threshold bin i at redshift z, and Poisson(/i, x) is the PDF of a Poisson process, x, of mean \x. 
A joint analysis of all critical sets, {£ k (v, z)}fc, is achieved via the product on k (assuming somewhat optimistically that all 
estimators are independent). 

We shall not carry such a maximization in full here. Instead, Figure 10 displays (very naive) likelihood contours for an 
experiment corresponding to the CFHT-LS extrapolation (Durret et al. 2011) of the cluster counts to a 1/4 sky experiment 
(2% error on D(a) for 40 bins, which corresponds to the appropriate rescaling of equation (88) given the prefactors of Table 4). 
Here we simply assume that the fit to the critical sets yields access to a noisy D(a\wo,w a ) via equation (89) (for a Gaussian 
noise with a SNR of 50= 1/0.02); for a grid of models obeying equation (89) we therefore compute the likelihood of a given 
draw. The corresponding 3-sigma contours suggest that critical sets could constrain Wo to ~ 5% and w a to ~ 10%. 

Beyond the obvious over-simplifications of this model, here are a few limitations/caveats which need to be discussed, (i) 
the accuracy of the predictions, £ (y, z) is as good as PT; it therefore requires a galaxy sample large enough to probe the larger 
scales, so that a 0.2. On scales of the order of 30 Mpc//i, given that cr(R, z) = ao(R)D(z), with ctq(R) the r.m.s of the cosmic 
density field at redshift zero smoothed over R, if we rely on the estimates for a given by equation (88), we may expect to probe 
D(z) near z ~ 0.5 13 ; (ii) for scale invariant power spectra, the cumulants entering equations (71), (76) and (81) and the shape 
parameter, 7 depend on the power-law index n s (see Appendix F2). For realistic ACDM cosmologies, following Bernardeau 
et al. (2002) one strategy could be to let the index n s be a (redshift dependent) parameter to be also estimated while fitting the 
counts 14 ; (iii) one could also estimate the optimal weighting strategy to best constrain (wo,w a ) using combinations of critical 
set estimators, varying threshold and redshifts weights. This is achieved by replacing Poisson (£^(z^, z), E^ A in equation (90) 
by Poisson {w z i k £\(yi, z), wf^E^A and minimizing the error on the estimated dark energy parameters w.r.t. the positive 
weights w* h ; (iv) in this paper, we ignored all issues regarding completeness, masking, etc. (see Section 5). In closing, note 
that the 2D statistics, equations (37), (47) and (63) could also be implemented on non-linear projected data, given some 
dynamical model for the projected cumulants, e.g. in the context of weak lensing (see e.g. Matsubara 2003; Yang et al. 2011; 
Kratochvil et al. 2011; Marian et al. 2011). 

Carrying out the road map sketched in this section while addressing these issues could be one of the target of the upcoming 
surveys that have been planned specifically to probe dark energy, either from ground-based facilities (eg VST-KIDS, DES, 
Pan- STARRS, LSST 15 ) or space-based observatories (EUCLID, SNAP and JDEM 16 ). 

In this paper, we have chosen to focus mainly on the non-Gaussianities arising from gravitational instability. It should be 
noted however that our formalism is more general. The Gram-Charlier expansion is valid for any mildly non-Gaussian field, 
while perturbation theory only enters the prediction of the cumulants. An important example of situation where the cumu- 
lants can be computed independently is the /nl parametrisation of non-Gaussianities in the CMB. In this parametrisation, 
cumulants can be computed and our formalism can be adapted to give predictions of the evolution of the properties of the 
field with the value of /nl- For example, Fedeli et al. (2011) recently used the skeleton to constrain /nl on two-dimensional 
weak-lensing maps. Their approach was numerical and our formalism could provide the theoretical framework for this type 
of investigation. 



5 CONCLUSION AND DISCUSSION 

We have introduced the Gram-Charlier expansion to describe the point probability distribution of a field and its derivatives 
(equations (11) and (24)) in the mildly non-Gaussian regime in two and three dimensions. We have derived from these 2D and 
3D PDF the non-Gaussian expressions for the Euler characteristic, equations (37) and (71), for differential extrema counts, 
equations (47) and (76), and for the skeleton differential length, equations (63) and (81). For the Minkowski functionals, 
including the Euler characteristic (and therefore the rare event tail of the differential number counts) the expansion is valid 
to arbitrary order in Gram-Charlier and is also given explicitly to third order in r.m.s of the density field. All differential 
distributions were checked against Monte Carlo realizations of the field and its critical sets: excellent agreement was found 

11 Here we assume for simplicity that the different counts are uncorrelated; an improvement would be to use multinomial statistics. 

12 As the Euler characteristic is cumulative, the model needs to be differentiated since equation (90) assumes that the counts are 
independent. In other words, £^ uicT (u,z) = d\jdv given by equation (48). 

13 Note that the invariants moments q 2 , Ji or J2 involve higher order derivatives, and are therefore sensitive to smaller, less linear scales. 

14 Alternatively, cumulants like equation (Fll) can be computed numerically for realistic power spectra. 

15 http: //www. astro-wise . org/projects/KIDS/, https : / /www. darkenergy survey . org/, http : //www. lsst . org/ 

16 http: //sci . esa. int/euclid/, http : //snap. lbl . gov/, http : //jdem. lbl .gov/ 
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a w'o - 1 



Figure 10. Left; D(a) versus a for the pseudo data together with a bundle of models with different values of wq and w a ; right: one, 
two, three, four and five sigmas likelihood contour in the (wo,w a ) plane for a one quarter sky experiment (consistent with a significant 
fraction of the number of clusters found in the CFHT-LS wide survey) leading to a 2 % Gaussian relative error on D(a) for 40 bins in a 
between 0.1 and 1. The red and blue dots correspond to the recovered most likely solution and the input value. For such an experiment, 
5 % on wo and 10 % on w a seems achievable. 



for extrema counts (and therefore for the Euler characteristic) for dispersion up to ~ 0.2. For the skeleton, some level of 
discrepancies still arises for some intermediate contrasts (see Figures 3 and 7). The origin of this disagreement may originate 
either from the stiff approximation behind the theory, or residual biases in the measurements. As argued in Section 4.2.4, the 
range of contrasts corresponding to this mismatch may in practice be masked out when fitting observed differential lengths. 
Furthermore, using gravitational clustering perturbation theory, we have shown how to predict all the 3D cumulants that 
appear in the Gram-Charlier expansion from non-linear gravitational dynamics. We can therefore predict the evolution of 
quantities such as the Euler characteristic, the density of extrema and the differential length of the skeleton as a function of 
r.m.s of the cosmic density field. The evolution of these critical sets can thus be related to the evolution of the underlying dark 
matter field. This connection is achieved via the gravitational distortion of quantities which can be measured independently 
of any monotonic bias. We have sketched a possible dark energy experiment based on the statistical analysis of these critical 
sets in order to constrain its equation of state, which in the very idealized regime we investigated, seems promising. 

Possible extension of this work beyond the scope of this paper involve i) implementing the dark energy estimation of 
section 4.2.4 first on realistic mock catalogs extracted from large scale simulations such as Teyssier et al. (2008), and eventually 
on the above cited surveys, while accounting for incomplete, magnitude limited samples and modeling the corresponding biases; 
one remaining task would be to demonstrate that the statistical analysis of geometrical critical sets can be competitive in 
practice to constrain this equation of state, compared to e.g. weak lensing (Kaiser 1998; Bridle & King 2007; Refregier et al. 
20f0; Pichon et al. 2009) or Baryonic acoustic oscillations probes (Zunckel, Gott & Lunnan 2011); ii) deriving the invariant 
JPDF for higher order derivatives and/or anti-derivatives (e.g. to link geometry to the kinematic of the flow along filaments), 
or for fields of higher dimensions (e.g. to study the saddle points of high dimensions landscapes), following the tracks presented 
in Pogosyan et al. (2009); iii) extending the theory of critical sets to walls and manifolds of higher dimensions; iv) exploring 
the implication of the departure from a Gaussian JPDF on other statistics, such as the mean connectivity of the peaks of 
the cosmic web, following Pichon et al. (2010); v) implementing second order PT, or the so-called renormalized perturbation 
theory (Crocce & Scoccimarro 2006) in order to extend the domain of applicability of these asymptotic expansions, and 
possibly build a generalized stable clustering model for our cumulants. vi) extending the formalism to N-points statistics, e.g. 
to investigate the non-Gaussian correlation of peaks, or non-Gaussian count-in-cells. 
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APPENDIX A: DERIVATIONS OF THE EULER CHARACTERISTIC IN 2 AND 3D 
Al Euler characteristic in 2D 

For a two-dimensional random field, the Euler characteristic xW) of the excursion set x > v is given by: 




(Al) 
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where P ex t is the distribution function under the condition of zero gradient: 



Pext(C,Ji, Ji) = 



dq 2 P((, q 2 , Ji, J2)Sd(xi)5d(x2). 



Using polar coordinates P ex t(C, Ji, J2) = ~P(G, Q 2 = 0, Ji, J2), and using the fact that Lj(0) = 1 for any j we have: 



P«xt(C,Ji,J2) = -GaD(C,Ji, Ja) 



00 i+2j+fc+2i=n 

i + E E 
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Remarkably, the Euler characteristic can be evaluated completely for this general distribution function. For this we note that 
I2 — \{J\ ~ J2) = \ (H2(Ji) + Li(J2)), thus, the calculations are reduced to the integration of products of Hermite and 
Laguerre polynomials. As the first step, the integral over J2 leaves only the terms with Lq(J2) (/ = 0) or Li(J2) (I = 1) in 
the Gram-Charlier expansion: 
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Let us first establish the important relation between some coefficients of the expansion for an isotropic field by evaluating 
x{— 00). Integrating over full range of £ selects i = terms, 
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and subsequent evaluation retains only j = n/2 — 1 in each even order of the expansion 
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The limiting Euler number \(— oo) is determined by the topology of the manifold which the sum of the cumulants in Eqs. (A6) 
must obey. We note that in this derivation we have obtained only the non-Gaussian contribution to x(~°°) which comes via 
hierarchy of the terms (q 23 12 ) , j > 1. The Gaussian term of the same structure corresponds to j — and was set to zero 

\ / GC 

by our normalization (^2) = (Ji) — (J2) = 0. It is an open question whether this is strictly consistent only with \{~ 00 ) — 
as, e.g., in an infinite Eucledian space, but requires more accurate treatment on the 2D sphere where \(~ 00 ) — 2. For now, 
we will leave the boundary term \(— 00) included in our expressions but unspecified. This is also justified by the practical 
case of the boundary topology introduced by a mask on the data. In this case the moments of the field remain the same as 
on unmasked manifold, but \(~ 00 ) measured in realizations will nevertherless reflect the topology of the mask, being now 
disentangled from moment statistics. 

Returning to the differential case, we use the following property of the Hermite polynomials: 
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Let us first treat the i = case, expanding an appearing product of Hermite polynomials back into Hermite polynomials 

Hp(Ji)Hk(Ji) = 7 — T-, rr-r, rrHk+ P -2s (Ji), the contribution to the Euler characteristic reads 
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We now use 



to obtain 
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where in the latter summations only terms that lead to non-negative index of Hermite polynomials should be included. The 
remaining i > terms of the Euler number calculation are 
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We need the following convolution property of the Hermite polynomials 



^ + 7J1 



7 



-fffc (Ji 



(A10) 



dJiG(Ji)exp - 



1 

to obtain 



xO)°>o = ^ — as=^ ex P 



(y + 7J1) 2 
2(1 -7 2 ) 



Hi - 



V + jJl \ rr I j V 1 / ^ 



(_!)«(! - 7 2 ) i 4 1 7*jr i+fc ( I ,) (All) 



00 i+2j+h=n r-^y+k 

E E ifTi 

n=3 i = l,j,k ' J ' 

00 i + 2j + fe + 2=n 

+ E E 

n = 3 i=l,j,k 



C g Ji ) (1-7 

GC 



2\V2 



min(2,fc) 

E 



2!-7 fe+2 ~ 2s 



s!(2 -*)!(&-«)! 
Cq 2J Ji k J2) (l - 7 2 ) l/2 ^Hi+k-i {v) 



Hi+k+i-2s (y) 



3=0 

i! j! it! 



(A12) 



We see that it has the same form as the second part of the Xi=o contribution which can be joined just by extending the range 
of i to start from zero. Thus, together with the result for \{— 00) we find the full Gram-Charlier expansion for the Euler 
characteristic, equation (37). 



A2 Euler characteristic in 3D 

Following the 2D calculation, one just averages —Is = — (Ji 3 — ZJ1J2 + 2J3) /27 over the full range of the second derivatives 

,3/2 



X{v) = — B3 / dJ a / _ dC / dJ 2 / _ , dJ 3 Pext(C, Ji, Ja, (A13) 



In 3D, we have 



o \3/2 

P«t(C,Ji,Ja,J8)= ( P(C,Ji,J 2 ,J 3 ). (A14) 



Writing ^3 = ^i?3( Ji) + |//i(Ji)Lj 3/ ^ 2 ' (| J2) + 2Jaj /27, one can perform all the integration explicitly using the orthogonality 
properties of the polynomial expansion of P ex t given by equation (24). First, integration over J 3 and J2 gives the intermediate 
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formula 
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where we have used l/ 1/2) (0) = (2j + l)!!/(2 J j!). The Euler number at the threshold that encompass all the manifold 
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As a function of the threshold, v, the integration is performed using the properties of the Hermite polynomials, first to 
get (given equations A7, A8, All) 
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and 



X(y)i>o = - wo3 , a / dJ 1 G(Ji)exp' 



27i?2 4tt 2 v-/— *- ^ 2(l- 7 2 ) 

00 i+2j' + A;— n 

E E ;> ,! /,i 

n—3 i=l,j,k 
00 -i + 2j + fc+2=n 



E e ^(-§y(cvww 



+ E E 

n—3 i=X,j,k 



z y (Cq 23 Ji k J2) H, 



i\ j\ fe! v 21 



v + 7J1 

^ + 7J1 



7J fc (Ji) 7J 3 (Ji) 



^(Ji)Hi(Ji) 



i + 2j + fc + 3=n 



+ E E 

n — 3 i— l.j.fc 



f + 7 Ji 



= (Jl) 



(A18) 



and, finally, equation (71). 



APPENDIX B: FITS TO EXTREMA COUNTS AND SKELETON IN 3D 
Bl Differential extrema counts 

The Gaussian distribution of maxima is a classical result and fitting formulas have been proposed e.g. in Bardeen et al. (1986). 
For the sake of consistency, we propose a formula using the same functional form as the next fits. For a scale-invariant power 
spectrum of spectral index n s , the Gaussian distribution of maxima is accurately fitted by : 
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The first order non Gaussian correction is fitted by: 
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For the saddle points, we have: 
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Replacing the field by its opposite (i.e. replacing v with —v and the odd Gram-Charlier coefficients by their opposite) yields 
the expressions for minima and the second kind of saddle points: 
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The corresponding fits for the first order corrections are thus obtained by changing the signs of the even Hermite polynomials. 
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B2 Skeleton length 

For n s £ [—2.5, 0], the Gaussian prediction for the differential skeleton length is well fitted by 
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APPENDIX C: COMPUTING THE SKELETON LENGTH IN THE HESSIAN EIGENFRAME 

The calculations made in Section 3.2.4 used the assumption that the properties of the components of the gradients are the 
same in the Hessian eigenframe as in any frame. This assumption translates into {x\) = {x\ 2 ) = 1/2, (i 2 ) = {x 2 2 } = 1/2, 
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Figure CI. Prediction for the 2D skeleton length using the components of the gradient in any frame (dashed line) or in the Hessian 
eigenframe (solid line). A small improvement is seen, especially near the maximum of the skeleton curve. However the correction remains 
very small, and its implementation also requires diagonalizing the Hessian and measuring two cumulants that will not be predicted by 
perturbation theory. 



and (Cil) = (C 2 ^ 2 ), {x2 2 Ji) = (x2 2 Ji)- These equalities are true in the Gaussian case: the first and second derivatives are 
independent so that the statistical properties of the gradient do not depend on the Hessian eigendirections. However it is no 
longer true when non-Gaussianities give rise to correlations between the gradient and the Hessian. To avoid this assumption, 
one can measure the values of (Cje!) and (x~2 2 Ji)- To do so the Hessian needs to be diagonalized at each point, which makes the 
measurements a bit more difficult. Furthermore the cumulants measured in the Hessian eigenframe seem difficult to predict 
via the perturbation theory described in the next section. 

The assumption on the variances (x\) and is more subtle. These terms are not explicitly present in the expansion. But 
they are actually implicitly assumed when the Gaussian kernel is made to built the Gram-Charlier expansion. If these variances 
are taken in a fixed frame, the isotropy and the normalization of the gradient (q 2 ) = 1 fix these variances to 1/2. But if the 
Hessian eigenframe is used, these variables are still Gaussian but their variance can change. The non-Gaussianities then affect 
the Gaussian kernel and the expression of the expansion. However this problem can be easily solved by using these variances 
to normalize the variables. We can use the variable Q — 1/ y/2{xl£)X2 instead of x^ to make sure that the variance is fixed and 
that the expansion using this new variable will not change. To summarize, if the properties of the gradient are not assumed 
to be the same in the Hessian eigenframe as in a fixed frame, one needs to diagonalize the Hessian, compute the components 
of the gradient in the eigenframe, normalize them to a fixed variance, measure the cumulants (CQ 2 ) and (Q 2 Ji) and use the 
previous analytical predictions with the new normalized variable Q instead of the original Xz. Figure CI shows the results of 
such a procedure. The difference with the previous results is almost negligible: the statistical properties of the components of 
the gradient are almost the same in the Hessian eigenframe as in a fixed frame and considering the apparition of a difference 
with non-Gaussianities is just a higher-order correction. 



APPENDIX D: CONNECTING THE EDGEWORTH TO THE GRAM-CHARLIER EXPANSION 

In this section, we give a formal derivation of the Edgeworth series (Blinnikov & Moessner 1998). For convenience, we restrict 
it to the case of a monovariate distribution, but it can be generalized directly (Kendall & Stuart 1977). Understanding the 
link between the Gram-Charlier and the Edgeworth expansions will be the key to a physically motivated resummation of the 
expansion. For a given PDF P(x), we define the moment generating function 



and the cumulant generating function 



/ + oo °° /„p\ 

P(x)e tx dx = ^j-t V , (Dl) 

°° It p \ 

C{t ) = E -^T tP - ( D2 ) 

p=2 P ' 

By definition of the cumulants, we have 

M(i) = exp (C(t)) . (D3) 
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Now, let's take the Fourier transform of the PDF: 



P(t) = — L / P(x)e" x dx. (D4) 



'2ty 
We have 

P(t) = M(it) = exp(C(it)). (D5) 
We can get the PDF by taking the inverse Fourier transform: 

P{x) = -}= / e" ite exp £ dt. (D6) 



In order to have an expansion around the Gaussian case, we extract the p = 2 term which corresponds to the purely Gaussian 
contribution: 

P( x ) = J= j +X e-'^'T' 2 exp ^p^ P J dt (D7) 
where a 2 = (a; 2 ). If we expand the last exponential, we will thus have the Gaussian term 
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The expansion can therefore be formally written as 
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Let us expand the exponential in Taylor series: 
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where G(aj) = exp(— x 2 /2a 2 ) /\ / '2na 2 is the Gaussian PDF. Developing powers of the operator allows us to rewrite equa- 
tion (Dll) as 
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By definition of the Hermite polynomials: (— 1) d G/dx = H k (x/a)G(x)/a it follows that 
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Joining the terms having the same Hermite polynomial would give the Gram-Char lier expansion. The Gram-Charlier coefficient 
{x n )cc can therefore be expressed in terms of cumulants: 

, Pl+P2=« ( P1) , p 2) 1 Pl+P2+P3=1 (x P1 ) ( X P2 ) (X P3 ) 

2! z — ' , Pi!p2! 3! ^ Pi!p2!p3! 

P1,P2=3 P1,P2,P3=3 

Pl+P2=™ 

where the notation means a summation on p\ and pi from 3 to 00 under the condition p\ + p2 = n. In the context 

Pl,J>2=3 

of gravitational instability, we have (x n ) oc a 2n ~ 2 . The expansion can be reordered according to the order of the cumulants, 
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introducing the coefficients S p = {x p ) c /a 



2p-2. 



p=3 / \Pl— 3p2— 3 
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The first terms are: 



- (i + + + + •••) + 

Ordering the terms w.r.t their order in a gives the Edgeworth expansion: 

P(x) = (l + a^Hzi-) + a 2 ( ^H 4 (-) + ^^H B (-)\ + a 3 ( ^H 5 (~) + ^ff 7 (-) + &L-H a (~)\ + .) G(x), 
K> V 6 ct V24 a' 72 a J J \120 V 144 a' 1296 K a'J J K ' 

which corresponds to the expression found in the literature (Bernardeau et al. 2002). Equation (D14) allows us to determine 
the orders to which a given coefficient contributes to the summation. Indeed, a coefficient with n fields can be decomposed as 
a sum of: 

• a linear (in cumulants) term, which is of order <j n ~ 2 , 

• if n > 3, a quadratic term, of order a pi ~ 2 + a P2 ~ 2 — a n ~ 4 , 

• if n > 3fc, a product of k cumulants, of order a n - 2fe - 2 
Thus the orders in a can be obtained as follow: 

• the a order comes from the n — 3 terms, 

• the a 2 order comes from the n = 4 term and part of the n — 6 one, 

• the a 3 order comes from the n = 5 term and part of the n — 7 and n — 9 ones, 

• the a 4 order comes from part of the n = 6, n = 8, n = 10 and n = 12 terms, 

• ... 

In the next section, we will see how to compute explicitly the required elements of the expansion up to some order in a. 



APPENDIX E: CRITICAL SETS IN THE EDGEWORTH PT EXPANSION 
El Truncating the Gram-Charlier expansion to order a 2 

According to Appendix D, the a 2 term arises from the 4th order cumulants and from part of the 6th order Gram-Charlier 
coefficients. Indeed, the 6th order coefficients contain quadratic S^-like terms. For example, Equation D14 gives (a; 6 } G c = 
(x 6 )c + 10(x 3 ) 2 . In general, these contributions can also be computed by expressing the Gram-Charlier coefficients in terms 
of moments using the expressions of the orthogonal polynomials (see equations (13), (25) and (26)), and then expressing 
moments in terms of cumulants. Perturbation theory ensures us that a term involving a cumulant with p fields will be of order 
a p ~ 217 . For example, the 6th order Hermite polynomial is : fle(C) = C 6 - 15C 4 + 45C 2 - 15. So: 

(C 6 )gc = <#6(C)> = <C 6 >-15<C 4 > + 30, 
where we have used the fact that (£ 2 ) = 1. We can then use the decomposition of moments into cumulants: 

<C 4 > = <C 4 ) C + 3<C 2 ) 2 , <C 6 > = <C 6 )c + 15<C 4 MC 2 > + io<C 3 ) 2 + 15<C 2 ) 3 - 

In the end, 

(C 6 )gc = <C 6 ) c + io(c 3 ) 2 . 

The first term is of order a 4 and the second one is of order a 2 . Table El give all a 2 terms coming from the 6th order 
Gram-Charlier coefficients used in the computation of the Euler characteristic. The expansion to 9th order is available online 
(http : //www. iap . f r /users/pi chon/Gram). 



17 Or a 2p 2 if not normalized. 

© 0000 RAS, MNRAS 000, 000-000 



34 C. Gay, C. Pichon and D. Pogosyan 



Coefficient 


a 2 


<C b )gc 


10<C S > 2 


<C b ^l>GC 


io<c :i >(C 2 Ji> 


<C 4 ^>GC 


6(C 2 JiX + 4<C 3 XC/i 2 } 


(C 3 ^1 2 >GC 


9<C 2 ^iXC^ 2 } + <C 3 XV> 


(C^lXcC 


6(CJi 2 > 2 +4(J 1 3 >(C 2 Ji> 


<CJl & >GC 


io(Ji a >(CJi 2 > 




10(,/i a } 2 


(CV)ao 


4<C 3 XC<X> 


(CVJl>oc 


(C 3 >to 2 ^i> + 3<C<? 2 XC 2 ./i> 


<CVJi 2 >gc 


2<Cq 2 Xai 2 > + 2(g 2 Ji><C 2 Ji> 


<C?^1 3 >GC 


<JXXCg 2 > + 3^JiXCJX> 


<9 2 Ji 4 >gc 


4<VX? 2 Ji> 


(g B >Gc 






Coefficient 




(C 4 ^2>GC 




(C^1^2>GC 


<C 3 XJiJ2> + 3<cj 2 XC 2 Ji> 


<C 2 ^1 2 ^2>GC 


2{C,J 2 ){C,Ji 2 ) + 2<JiJ 2 >(C 2 Ji> 


(C^2>GC 


<Ji 3 XC-/2> + 3<JiJ 2 >(C^i 2 > 


<Jl 4 J 2 >GC 


4<Ji 3 >(J 1 J 2 > 




2<C9 2 XC^2> 


(C^l^GC 


<C</ 2 >< u 'i^2> + <g 2 JiXC^2> 


/ 2 T 2 T \ 

(g 2 Jl J2>GC 


2<g 2 Ji)<JiJ2> 




(C 3 >(J 3 > 


(CVlJ 3 >GC 


<C 2 Jl><J3> 


(C^1 2 ^3>GC 


<ai 2 X^3> 


<Jl 3 J 3 >GC 


(Jl">(Js> 


(C<? 2 -/3>GC 


<C9 2 XJ3> 


(9 2 Jl^3>GC 





2D 



Coefficient 


a 2 


<cv>gc 


4<C9 2 X 


<C9 4 ^i>gc 


4(C<? 2 ><9 2 Ji> 


(9 4 Jl 2 >GC 


4<9 2 Ji> 2 



3D 



Coefficient 


a 2 


<fV>GC 


WX 


<C<? 4 ^1>GC 




<<M>GC 




(<? 4 ^2>GC 





Table El. a 2 order terms from the n = 6 Gram Charlier order. Only terms contributing to the Euler characteristic are shown. Terms 
involving J3 do not exist in 2D, and those involving q 4 differ in 2 and 3 dimensions. 



E2 Minkowski functionals 

E2. 1 2D length of isocontours 

In order to illustrate how to use Table El, let us consider the reordering of the length of the isocontours in 2 dimensions 
presented in Section 3.2.1. To obtain an expression up to a 2 order, one has to compute the n = 3, 4 and 6 term from the 
general expression (33): 



£W = 2XT 



1 + \{x Z )H 3 (v) + lix^Wu) + ±(x%m{ V ) + \(x 2 q 2 ) c H 2 {») - ±{ q %H (v) 

+ ^{x 6 )gcH 6 (u) + l(xV) GC H 4 M - ^L{ X 2 q% c H 2 (v) + ±{q 6 ) GC H (v)) • (El) 



Using the results from table El , this expression can be truncated at the a order: 

Ciy) = ^f"^ (l + \(z 3 )Hz{v) + -(zgViM + ^(x A )Mv) + \{x 2 q 2 ) c H 2 {v) - ~{q\H a {v) 

+ ±{x 3 ) 2 H 6 (v) + ±(x 3 ){xq 2 )H 4 (v) - ±-{xq 2 ) 2 H 2 {v)\. (E2) 



72 



12 



32 



E2.2 3D surface of the isocontours 

Following the same route, equation (68) in Section 3.3.1 becomes to the required order 



2 

\/3V 



A(y) = -ire-^ I 1 + hx^)H!(u) + \(x s )H,{u) - ^) c H (u) + hx 2 q 2 ) c H 2 {v) + l(x%H*{u) 



+ 



40 
9 



560 



27 (q 6 )ocH (u)^^(x 2 q 4 ) GC H 2 (u) + ^(x 4 q 2 )ocH 4 (u)+ ' ' 



240 



) GC H s (v) 



(E3) 



Since no second derivative is involved in this simple example, x can be used instead of £. 
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Using the expression of the a contribution of the n = 6 terms, we get the a expression, and its truncation to a order reads: 



1 + \(xq % )H x {v) + \(x*)H 3 {v) - ^(, 4 ) c ffoH + \(x 2 q 2 ) c H 2 {v) + ~(x A ) c H 4 (u) 



(xq 2 } 2 H 2 (u) + \{xq 2 ){x 3 )H,{v) + ±(x 3 ) 2 H a (u) j . (E U 



All other terms in equation (E3) contribute to orders in a higher than 2. 



E3 2D Euler characteristic 

It is interesting to look at the first correction to the Euler characteristic since it has been predicted in other papers using a 
somewhat different approach (Matsubara 1994, 2003). The a order term is given by the n = 3 term of the Gram-Charlier 
expansion. In 2 dimensions, equation (37) yields 



4nV2nR'i 



7 2 HiH + (2j(q 2 Ji) + v/l-7 2 (C J i 2 > -7<J'i 3 ) - \/l-7 2 <C^)+7(.AJ 2 )) Ho(v) 
+ (-7 2 V^7 i (Cq 2 ) + 7 3 (q 2 Ji) - 7(1 - 7 2 ) <C 2 Ji) + 2j 2 ^l^f(CJ 1 2 ) - 7 3 (J! 3 )) H 2 {v) 

+ (y (i - 7 2 ) 3/2 <C 3 > + £(i - 7 2 ) <C 2 Jx> + ^x/W (C^ 2 ) - £ < Ji 3 )) 

This result can be simplified using the variable x and 7 2 : 
1 



7 2 ffiM + (2 7 (q 2 h) + 4 (a;Ja» flbM - {j 2 (xq 2 ) + 7 (z 2 /i» H*{v) + jr- (x 3 ) H 4 (u) 



4tvV2itR 2 

The expression given by Matsubara (2003) for the two-dimensional genus is 



e 2 



/C(0) no(l) o(2) 



where 



5 



(o) 



(P 3 ) 



,(D_ 3 (p 2 V 2 p) 

4 crVi 2 : 



7 (2) 



3d <(Vp) 2 V 2 p) 



(E5) 



(E6) 



(E7) 



(E8) 



a* 4 a^ai* 2{d — 1) cti 4 

We use here the unnormalized field, p, since Matsubara has a different normalization convention for the derivatives. With our 
normalization and for d — 2, these quantities can be expressed in terms of the rotational invariants: 



5 (0) = 



<*») 



3 (x 2 h) 



(E9) 



ct 4 70" 

Factorizing 7 2 from our expression for the Euler characteristic gives the same prefactor since y/R* = o\ja. The Hq(i>) term 
is obviously the same: 

5 (0) 1 , 3^ 



a = -(x 3 ). 

6 6 W 



Using the isotropy relations given in Matsubara (2003), one can check the identities: 



(xq 2 ) = (x 2 /!). 



{xh) = -~-y(q 2 Ii) : 



(E10) 



(Ell) 



which are needed to prove that the coefficients of H 2 (y) and H±{v) agree. Similarly one can express the Euler characteristic 
up to <7 2 order. The result is more compact when expressed in terms of x and I 2 : 

XH = ^J^ R f ~^ h'tfiH + (2 7 <g 2 /i> +4(^/2)) H {v) - ( 7 2 (xq 2 )+ 7 (x 2 h)) H 2 (v) + ^ (x 3 ) H 4 (u) 



+ 



-<g 4 ) c + 2 1 (xq 2 J 1 ) c + 2(x 2 I 2 )^j H^v) - (y <*¥}c + ^ 3 Ji>c) H a (u) + ^(x 4 ) c H 5 (u) 
- (4 7 (xg 2 )(g 2 J 1 ) + (q 2 J 1 )(x 2 J 1 ) + 4(xq 2 )(xl 2 }) H x {v) 
+ (^(xq 2 ) 2 + l(x 3 )(q 2 .h) + 1 (xq 2 )(x 2 J 1 ) + \(x 2 -h) 2 + l(x 3 ){xl 2 )) H 3 {v) 



(l(xq 2 }(x 3 ) + J^XsVi)) H s (v) + ^(x 3 ) 2 ^^)). (E12) 
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This result can be simplified using the isotropy relations (Ell): 



xiy) 



- I y(9 )c + 27(i9 Ji)c + 2{x 2 I 2 ) c ) H x {v) 



^-(x 2 q 2 ) c + I^JrJc) ffsW + ^(^)cH 5 H 



12 



\{q 2 J,){x 2 J^{v) - l{x*){ q 2 J,)H 3 {v) - ^(x 3 )(x 2 Ji)H 6 {v) + ^(x 3 ) 2 H,i,n 



This expression agrees with Matsubara (2010) given the following identities: 



1 



(x 2 q 2 ) c = -^(x 3 J 1 ) c 
37 



(E13) 



(E14) 



<g 4 ) c + -{x q 2 Jx) c + ^Vl 2 >c = . 
7 7 

These relations are consistent with the measurements. The compactness of the expansion in invariants allows us to compute 
the a 3 correction: 



X (3) H = (ll(q 4 )c(q 2 Jl) ~7(<M>c + ^(q 2 J 1 ) 2 (x 2 J 1 ) - ^(q 2 h}(x 2 J 1 ) + ~ 7 <g 4 ) c (q 4 J,) + 4(q 2 {xq 2 Jx) c - Mxq 2 I 2 ) c ) H (u) 
\2 772 J 



-i- ( 5<x 3 )( 8 a Ji) a - -(q'JiHx'JJc - ^(q'J^qX + ^{q%(x* Jr) + y(^Jt) c - -{q'J^Jtf + (x 2 J ,) (xq 2 J,) 



7/2 



2 2\ ,7/4 



+ i<:r 2 /2>c(2- 2 J 1 > - ^(g 2 /2)c<^ 2 ) + ^(x 9 4 )c - ^7 2 <9 4 )c(^ 4 )c + §<* 3 Ja>c) fl- a (f)+ 



+ ( "T^W^c - ^7(^>c<g 2 Ji) - ~(x 3 )(q 2 Jl}(x 2 Ji} - ^7(^Ji)c + ^1 2 {q 4 )c{x 3 } - + l{x 3 ){x 2 I 2 )c ) 



12 



12 



12 



- I 7^7 V>c - ^l{x 3 ) 2 {q 2 A) - ^l{x 3 ){x 3 A) c - l 7 (x 4 ) c ( 2; 2 J 1 ) - i 7 V>(z¥)c ) //,.„., 



120 



48 
1 

144 



+ 1 m j2{x3){x4)c ~ llJ^ 3 ) 2 ^ 2J l))^( ^/ ) + i46 72<a:2>3Jffl ' ll '' , ,Kr> ' 



E4 3D Euler characteristic 

In 3D, the first-order term from equation (71) is 
1 



27 Rl \2n 



3\ 3/2 1 



2vr 



+ 



+ ^7 2 (q 2 Ji) + 3 7 \/r^<CJi 2 ) - 3 7 2 < Ji 3 ) - 3 7 \/r 3 7 ? (C^) + 3 7 2 <JiJ 2 >) 
(-^a/W <C? 2 ) + §7* <a 2 Ji> - ^7 2 (1 ~7 2 ) <C 2 Ji> +37 3 \^^<CJi 2 > ~ §7* (Ji 3 )) H 3 (v) 
+ (^7 3 (1 " 7 2 ) 3/2 <C 3 ) ~ ^7 4 (1 ~ 7 2 ) <C 2 Ji) + \l 5 V^T 2 (CJi 2 ) - \^ (J! 3 )) H«M 



As in 2D, it can be simplified by using the variables x, I 2 and I3 



xW) = 



1 



27 Rl \2tv 



3 \ 3/2 1 



2tt 



J S H 2 (v) + ( ^7 2 (q 2 h) + 9 7 (xh) ) fliW 



(§7 3 (xq 2 ) + l-y 2 <* 2 /i>) H 3 (u) + i 7 3 <x 3 ) HsM 



Matsubara's expression for the three dimensional genus is 

3 



Gs(f) 



1 / 0-1 
(2tt)2 VV3a 



/5(D) 



H 2 (v) + l^~H s (y) + S w H 3 (u) + S^H^u) ) a 



with 



S m = - (x 3 ) , 



g(1)= 3 (s 2 /i) _ 
4 7<r 



5 (2) = 



9 (q 2 h) 



(E16) 



(E17) 



(E18) 



(E19) 



a ' 4 70" 4 70" 

Factorizing the expression for %( v ) with 7 3 leads to the same prefactor except for the sign. This sign difference comes from the 
fact that xiy) ls the Euler characteristic of the excursion set while Gz(y) is the genus of the isocontour and thus x( v ) = ~Gz{v), 
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Using the isotropy identities given in equation (Ell), it is then easy to show that each term coincides. Finally, the 3D Euler 
characteristic up to a 2 order expressed in terms ol x, I2 and 73 is: 

+ \l Z (x Z )H,{v) - (~"f(q 2 I 2 )c + 27(x/ 3 )c) H (u) + (~7V>« + Ji>c + ~7<^2>c) H a (u) 

~ (|7 3 <^g 2 )c + ^tV^c) H*{u) + ^{x 4 ) c H 6 {u) + ^ 7 ( g 2 J 1 } 2 + ^-{q 2 J^xh) + y <^ 2 ><I 3 >) H (u) 
- (f 7W><<M> + ^7(g 2 Ji)(^Ji> + y 7 (^ 2 )(x/ 2 ) + \{x 2 Ji){xh) + ^ 3 / 3 >) H 2 {v) 
+ (j 7 3 <^ 2 ) 2 + | 7 V>(<M> + ^WX^i) + ^7(^Ji> 2 + ~7<* 3 >(z/2>) H 4 (u) 



1 3/ 2\ / 3\ , 1 ,2 / 3\ / 2 T \ A tt f \ 1 1 3/ 3\2 



rWIM + 7 7 2 (^><^Ji) ) Heiu) + ^( x yHsiv) 



(E20) 



APPENDIX F: GEOMETRIC CUMULANTS FROM PERTURBATION THEORY 

In this section, we will compute the 3rd order cumulants involved in the non-Gaussian expansion in the case where non- 
Gaussianities arise from the non-linear gravitational evolution. 

Fl Derivation for an EdS universe 

For a three dimensional field, the cumulants involving the field and its derivatives can be written as 

Ud^d^d^S) (d^d^d^S) [d^d^d^S)^ (Fl) 

where df is the n-th derivative with respect to the i-th coordinate. Using gravitational perturbation theory (Bernardeau et al. 
2002), we expand the field as S — S^- 1 ' + + . . . . Inserting this expansion in equation (Fl) gives the lowest non-vanishing 
order contribution to the cumulant: 

( (9f 1 d% 2 9g 3 5) (9f 1 d? 2 2 d§ 3 8) (c\ 71 <9 2 72 dj 3 <*)) = ( (di 1 da 2 ^ 3 8 (1) ) (df 1 d§ 2 d§ 3 5 (1) ) {dj 1 d] 2 c^ 3 S <2) )^ 

+ /(a 1 ai a^ 3 a 3 l3 (5 (1) )(af i a| 2 af 3 <5 (2) )(a7 i aj 2 a3 y3 5 (1) )^ + ((d^d^d^s^xd^d^d^s^^id^d^d^s^)^ . (F2) 

Let us compute the first term in equation (F2), 1 = {{d^d™ 2 d^ 3 6 (1) )idf 1 d§ 2 d§ 3 6 i - 1) )idJ 1 d^d^5 (2) )), which is best expressed 
in Fourier space as (with k'- 7 ' the jth component of vector k) 

i=|d 3 k 1 d s k 2 d 3 k3 e i(ki+k2+k3) ' x ((zkW) Qi (iki 2, ) Q2 (i^yu^i^) 

^Y 1 (i^y* (i^y ^) (k2) ^iy ^y ( ik [ai)- 5 (2) (k3) ^ . (F3) 

In equation (F3), we assumed zero smoothing of the field for now. Using the notations of Bernardeau et al. (2002), given the 
Euler and continuity equation for an EdS universe, we can express <^ 2 ' in terms of tr 1 ' as : 

<5 (2) (k) = /d 3 k 1 d 3 k 2 «5 I5 (k-k 1 -k 2 )F 2 (k 1 ,k 2 )5 (1) (k 1 )<5 (1) (k 2 ), with f a (ki, k 2 ) = I + + % (kl a!" 2 / ■ ( F4 ) 

J < K\ ( fel K2 

Thus equation (F4) becomes: 

1= j d 3 k 1 d 3 k 2 d 3 k 4 d 3 k 5e ^ +k2+k -+ k ^- (zkW)" 1 (i^y 2 ( lk f)) Q3 (ikW)" 1 (zkl 21 )" 2 (z4 3] )" 3 

(i(kW+kW)) T1 (i((^ (F5) 
To simplify, we can combine the derivatives and F2 into a generalized shape factor: 

,/?3 



1)k2! k 4l k 5) = (* k w) ai (z^)" 2 (z^)" 3 ( ik w) ft (i^y (on)' 

(i (kW +k W)) 71 (< ((kf +4 2] )) 72 (< (kf +kl 3 ))) 73 F 2 (k 4 ,k 5 ), (F6) 

which gives 

1= y d 3 k 1 d 3 k a d 3 k 4 d 3 k 8 e i(kl+ka+ ^ +kB ^^ (F7) 
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The integrated moment is then expanded using Wick's theorem: 

( ( 5( 1 )(k 1 )5 (1) (k 2 )5( 1 )(k 4 )«5( 1 )(k 5 )) = (^(kO^Oca)) (5 (1) (k 4 )5«(k 5 )) + 

(^ 1 )(k 1 )5 (1) (k 4 )) (5 (1) (k 2 )<5«(k 5 )) + («5< 1 »(k 1 )5 (1) (k 5 )) (^ 1) (k 2 )^ 1 )(k 4 )) . (F8) 
The first term leads to P 2 (k 4 , — k 4 ) which vanishes, while the second and third terms are equivalent: 

X = 2^ d 3 k 1 d 3 k 2 d 3 k 4 d 3 k 5 J : Q , / 3, 7 (k 1 ,k 2 ,k 4 ,k5)(5 I5 (k 1 +k 4 )P(fci)5 D (k 2 +k 5 )P(fc 2 ). (F9) 
Integrating over k 4 and ks gives: 

X = 2^ d 3 k 1 d 3 k 2 ^^, 7 (k 1 ,k 2 )P(fc 1 )P(fc 2 ), (FIO) 

where we use the notation T a ,p,~i (ki, k 2 ) = Fa,p,~tQz-i, k 2 , — ki, — k 2 ) (see equation (86) in the main text). Finally, combining 
the three terms in equation (F2) gives us the cumulant: 

{{d^d^d^S) (d^d^dpS) (Bp8?a2>S)) = 2 j d 3 k!d 3 k 2 (^.fl^tk^k,) + ^.^(ki.kj) + JF 7 , a , /3 (k 1 ,k 2 ))P(fc 1 )P(fe 2 ). 

Let us now take into account the smoothing of the field over a scale, R; the cumulants become: 

2 j d 3 k!d 3 k 2 tF"a,/j, 7 (ki, k 2 ) + J>, 7 , Q (ki, ka) + TwjQlx, k 2 )] P(fc 1 )P(fc 2 )Vy(fciP)Vy(fc 2 P)Vy(|k 1 + k 2 |P) . (Fll) 

Equation (Fll) is the general expression for the cumulants used in the main text. The appearance of the dependence on the 
relative orientation of the wave vectors ki and k 2 in the filter factor W(|ki + k 2 |P) is the source of most of the complexity 
of the theory and, in some sense, the essence of perturbation theory. It reflects the fact that the nonlinear field, smoothed at 
radius R, is not determined solely by the average quantities at this radius, but also but what happens at shorter scales. For 
Gaussian filtering, the filter function can be expanded in Legendre series with respect to ki • k 2 (Lokas et al. 1994): 



W(\k 1+ k 2 \R) = cxp (-^kl + k 2 2 )R 2 j Y,(-l) e (2e+l)P e C^-f J /, +1 / 2 (fcifc2P 2 )^^^. (F12) 

Then, ^F a ,p,t can be decomposed on the basis of the Legendre polynomials and the integration over the angles then just 
requires the orthogonality relation of the Legendre polynomials: 

1 2 

dxP m {x)P n (x) = — -5 m ,n. (F13) 

! 2m + 1 

It also means that in practice, one does not need Legendre polynomials of order higher than the degree of the integrated term 
(here, 2 from P 2 plus the number of derivatives) and can truncate the expansion. The result is an integral of Bessel functions, 
which can be expressed using hypergeometric functions. The results for all the independent third order cumulants is given in 
Appendix F2. 



/_ 



F2 Gaussian Altered scale invariant geometric cumulants 

Using the above defined procedure, the cumulants can be computed for a Gaussian filter. For a scale-invariant power- 
spectrum of index n (called n s in the main text), the results can be analytically expressed using the hypergeometric function 
2 Pi(a, b, c, x). For example, the result for the skewness is already known (Lokas et al. 1994): 

1, 3 X o t? /3 + n 3 + n 3 l\ 1 . . /3 + n 3 + n 5 l\ . . 

-<* > = 3 afi [ — , — , 2 - 4 ) 7(8 + 7n) 2 Pi { — , 2' 4J ' ( F14 ) 

Proceeding the same way with the cumulants involving derivatives of the field yields: 

1. 2 , 4(48 + 62n + 2ln 2 ) ^ /3 + n 3 + n 3 l\ 6(3 + n)(8 + 7n) /3 + n 5 + n 3 l\ 

a {XXl > = 21* 2Fl { — ' — ' 2' 4 J 21^ 2Fl { — ' — ' 2 ' 4 J ' (F15) 

1. 2, 4(30720 + 51456n + 37092n 2 + 13828n 3 + 2579n 4 + 189n 5 ) ( fl „ , ^ / 3 + n 3 + n 11 
a {XXl1 } = 35n 2 (2 + n) 2 (4 + n) 2 (5 + n)3 |^ (4 + 2n) aFl { — > ~ ' ~2 ' 4 

,, /3 + n 5 + n 1 1\ \ 3(3840 + 5016n + 2748n 2 + 693n 3 + 63n 4 ) 

+ 3 2^1 ~ , ^ ,— ^: T + 



2'4J j 35n 2 (2 + n)(4 + n)(5 + n)(6 + n) 
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1. , 4(-46080 - 37984n - 2848n 2 + 5718n 3 + 1949n 4 + 189n 5 ) „, ^/3 + n 3 + n 1 V 

a <^ 22 > = 105n 2 (2 + n) 2 (4 + n) 2 (5 + n) 2 ^ + 2n) ^ {^T> ^T> ~ 2 ' 4, 

,, ' 3 + n 5 + n _1 l\ \ 3(-5760 - 2624n + 608n 2 + 476n 3 + 63n 4 ) 
+ 2 "' l| 2 ' 2 ' 2'4j) + 105n 2 (2 + n)(4 + n)(5 + n)(6 + n) 

, r , . 3 + n 7 + n 1 1\ /3 + n 7 + n 1 1 

3 a Fi ( 2 ' 2 ' _ 2 ' 4 / ^ n 2 1 I —'—'2' 4 



1. 2 , 4(9120 + 96176n + 57062n 2 + 17883n 3 + 2894n 4 + 189n 5 ) /.„ . ^ /3 + n 3 + n 11 
* {XX12 > = 105n 2 (2 + n) 2 (4 + n) 2 (5 + n) 2 I (4 + " } 2Fl ^'"2' 4 



+ 3 2 Fi 



3+n5+n 11 
2 ' 2 '~2' 4 



3(8640 + 8836n + 3818n 2 + 791n 3 + 63n 4 ) 
105n 2 (2 + n)(4 + n)(5 + n)(6 + n) 



, , 3 + n 7 + n 1 1\ , „ /3 + n 7 + n 1 1 
3 2 Fx ( 2 ' 2 ' _ 2 ' 4 J 12 l^-'-2~'2'4 



The last ones also depend on the spectral parameter 7 = \/(n + 3)/(n + 5): 



-{xiX 2 Xy 2 } 



42n 2 (2 + n) 



(96 + 50n + 7n 2 ) 2 F 1 '' 5 + n 5 + n 3 1 



;(a;i 2 ^i3^ 2 3) = 



87 



245n 2 (2 + n)(4 + n) \ ^° + ^ + ^ + 395 " S + 2 ^ ^ I' I 



(4320 



+ 3438n + 1180n 2 + 201n 3 + 14n 4 ) 2 F, (*±^, 1±±, |, 1 



47 



a< W > = ~ 735^(2 + n)(4 + n) ( 4(256 °° + 2128 °^ + 8 ° 02 " 2 + 1557nS + ™^ ^ ^ 2 ' 4 

- 3(19200 + 15280n + 5582n 2 + 1059n 3 + 54n 4 ) 2 F 1 ( |, j 



1, 2 > 47 

ff ^11x23 ; 735n 2 (2 + n)(4 + n) ^ 



8(1280 + 1064n + 608n 2 + 186n 3 + 21n 4 ) 2 Fi 
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;(£ii2 22 :r33} = - 
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+ (11520 + 9168n + 2134n 2 + 39n 3 - 28n 4 ) 2 Fi ^ 



2 ' 2 ' 2' 4 
5+n 7+n 3 1 



2 ' 2' 4 



{xii 2 x 22 ) 



47 

~735n 2 (2 + n)(4 + n) 



4(-1024- 8512n- 706n 2 + 675n 3 + 126n 4 ) 2 Fi ( 1 



- 3(7680 + 6112n + 410n 2 - 471n 3 - 84n 4 ) 2 Fi 



2 ' 2 ' 2' 4 
5+n 7+n 3 1 



2 ' 2' 4 



-(Xl X 22 ) = r 

<t 21n 2 (2 + n) 



CIS + 200, ^:,/:( :! 2 ".%". 3 .; 
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5 + n 5 + n 3 1 
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n s =0 


n s 


-1 


n s 


= -2 




prediction 


measurement 


prediction 


measurement 


prediction 


measurement 


(x°)/<J 


3.144 


3.08 ±0.09 


3.468 


3.53 ±0.15 


4.022 


4.3 ±0.4 


(xx-i 2 )/a 


0.699 


0.68 ±0.01 


0.771 


0.795 ±0.032 


0.894 


1.0 ±0.1 


(xxn 2 )/a 


0.567 


0.554 ± 0.009 


0.612 


0.636 ±0.026 


0.661 


0.8 ±0.1 


(xx 11 x 22 )/a 


0.361 


0.348 ± 0.006 


0.340 


0.361 ± 0.018 


0.311 


0.44 ± 0.06 


(xxi2 2 )/a 


0.103 


0.102 ±0.002 


0.136 


0.137 ± 0.005 


0.175 


0.20 ±0.02 


(xiX 2 x 12 )/(7 


0.111 


0.107 ± 0.002 


0.096 


0.106 ±0.007 


0.0784 


0.13 ±0.02 


{xi 2 xnx 2 z}/(T 


0.0124 


0.0122 ±0.0003 


0.007 


0.010 ±0.01 


0.00254 


0.009 ± 0.003 


(xnx 12 2 )/a 


-0.0068 


-0.0061 ± 0.0006 


-0.016 


-0.014 ±0.001 


-0.0237 


-0.025 ± 0.003 


(xiix 2 3 2 ) / 'a 


-0.0316 


-0.0306 ± 0.0006 


-0.031 


-0.034 ±0.002 


-0.0288 


-0.046 ±0.007 


(xnx 22 x 33 )/(7 


-0.120 


-0.116 ±0.002 


-0.108 


-0.121 ± 0.008 


-0.0914 


-0.15 ± 0.03 


(xii 2 x 22 )/a 


-0.183 


-0.177 ± 0.004 


-0.170 


-0.188 ±0.012 


-0.149 


-0.24 ±0.04 


{xi 2 x 22 )/a 


-0.222 


-0.214 ± 0.005 


-0.193 


-0.213 ±0.014 


-0.157 


-0.25 ± 0.04 


(xn :i )/a 


-0.210 


-0.203 ±0.004 


-0.235 


-0.244 ±0.011 


-0.244 


-0.35 ± 0.05 


{x 2 xii}/a 


-1.08 


-1.05 ±0.02 


-1.090 


-1.097 ±0.055 


-1.03 


-1.3 ±0.1 



Table Fl. Predicted and measured (Sec. 4.2.3.) third order 3D coordinate cumulants. According to perturbation theory, the dimensionless 
third order cumulants are proportional to a. 



- (23040 + 18336n ± 7306n 2 ± 1569n 3 ± 140n 4 ) 2 fi (^^>-^r^>| ) l> ( F26 ) 



All other cumulants are either identical to one of these because of isotropy, or null (mostly because of the parity of the 
derivatives). These formulae are thus sufficient to build any combination (e.g. the cumulants of the rotational invariants). In 
Section 4.2.3, Tables Fl and 3 give numerical values and show they agree with measurements. 



APPENDIX G: EULER CHARACTERISTIC ALGORITHM 

The very fast, linear in the number of pixels, Euler2D algorithm that we use to measure the Euler characteristic in 2D 
is based on the Gauss-Bonnet theorem. According to this theorem, the Euler characteristic of a region can be obtained by 
integrating the curvature over the boundary surface. We however note that explicit integration of the curvature is not needed, 
since the any continuous deformation of the regions conserve the topology. Hence, on a grid the curvature can be thought of 
as concentrated in the corners of the cells at the boundary between the cells above and below the threshold and just counted 
by addition. What is needed is to assign the appropriate curvature weight for each boundary grid vertex, depending on which 
cells that form it are below and which are above the threshold. The sum of the weights over all the vertices will give us the 
Euler characteristic of the excursion set. Our implementation of the algorithm is complemented by the cluster analysis which 
allow to count the Euler characteristics of the individual isolated regions. 

To determine the weights, we use a bootstrap method: we consider elementary configurations for which the Euler char- 
acteristic is known in order to determine the number to associate with a new geometrical configuration. 

For a 2D regular grid, the geometric configurations can be classified in six categories, represented in Figure G2. The 
number associated with the vertex in the first configuration (when all four surrounding pixels are not part of a region) is 
trivial: it is zero, since adding some empty pixels around a region does not change its topology. Similarly a vertex surrounded 
by four pixels has also the weight of zero since vertices interior to the excursion set does not change its Euler characteristic. 
The second configuration can be studied by looking at the case of a region made of only one isolated pixel (Figure Gl(a)): its 
Euler characteristic is unity and all the four corners of the pixel are equivalent and should thus carry the weight 1/4. These 
vertices are of a type when three neighbours are outside the region and one is inside. Recording this value, we can know look 
at the next category of pixel. Figure Gl(b) shows that if the region is made of two adjacent pixels, the Euler characteristic is 
still unity and the intermediate corners of pixels must contribute the weight of zero since four corner vertices are of the type 
considered before and carry the weight of 1/4. We can then consider a L-shaped region made of three pixels (Figure Gl(c)) 
to find that a vertex surrounded by three pixels in a region has a weight of —1/4. The last remaining configuration is more 
problematic: when two pixels touch only by a corner (Figure Gl(d)), one has to decide whether they are part of the same 
region or not. Depending on this choice, the weight is 1/2 or —1/2. In case we do track the topology of individual regions one 
of the decisions has to be made, which will bias the statistical result (but will be exact for a given definition of a connected 
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1/4 



1/4 1/4 1/4 1/4 



(a) X = 1 (b) X = 1 (c) X = 1 (d) X = 1 or X = 2 

Figure Gl. Simple configurations used to determine the number of the geometrical configurations. The Euler characteristic of the 
example is written below the drawings. The red numbers are the new weights that are determined at the current step of the bootstrap, 
while the blue ones are the previously determined ones. 

BB B B K B ■ 

i o -\ 

Figure G2. Examples of the 6 geometrical configurations that can appear with a 2D regular grid. The grey pixels are part of the region 
while the white ones are outside of it. The numbers indicate the curvature weight with "statistical" choice for connectivity in the case of 
touching vertices (the fourth plot). 

region). If we are interested in statistical analysis of the total Euler characteristic of the complete excursion set we note that, 
statistically in 2D, the two pixels touching by just the vertex are as likely to be connected as they are not. Thus, we make a 
statistical choice of assigning the weight zero to this last vertex type. 

The final assignment of curvature weights is given in the Figure G2. One can easily repeat the procedure in 3D, or 
generalize it onto more complicated grids. It can be noticed that for the simple case of a 2D regular grid, the result is simply 
the number of turns the boundary contour rotates at a given corner (Weinberg 1988). Furthermore, with our statistical 
prescription, the weights associated with different geometrical configurations only depends on the number of surrounding 
pixels in the region, not on their configuration, which makes coding the algorithm very simple. 

In complete implementation the genus algorithm is complimented by the cluster analysis, so that the Euler characteristic 
of each individual isolated cluster can be tracked. This approach has also been implemented for 3D grids and more complicated 
HEALPix pixelization in Euler3D and EulerHealpix families of codes. The codes are available from the authors. 
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